We perform here different aspects of beta diversity analysis of the microbiota data gathered from fish gill mucus and bacterioplankton. Most relevant sections are the distance based 6.3.1 PCOA and 6.4.3 RDA analyses as well as the 6.3.5 Mantel test for global matrix correlations. Biomarker discovery via multipratt gives some indications for taxa specific for fish, bacterioplankton and seasonal sampling. Overall it seems adviced to run several biomarker discovery techniques in parallel but in this variable dataset they give little insides and are more meant as a backup for the network analysis.

1 Global Setup

1.1 Path

before you start: .rs.restartR()

1.2 Packages

1.3 Functions

1.4 Input Files

1.5 Tutorials

1.6 Setup Analysis

#-

6 Beta Diversity

6.1 PCA

CLR tranformed data in PCA analysis

##################
# SampleDist_PCAs# Publication 
##################
#pslist <- pslistraw
for (i in names(pslist)[grepl("TSEc.l.r", names(pslist))]){
  require(plyr)
  require(ggrepel) 
  require(cowplot)
  require(DESeq2)
  require(matrixStats)
  TSE <- pslist[[i]]
  tryCatch({
  g       <- paste(i); print(i) 
  gg      <- sub('TSEc.l.r_', '', g)
  
  COL <- col.Palette$col.Palette.RepsSpecs[names(col.Palette$col.Palette.RepsSpecs) %in% TSE@colData$RepsSpecs]
  
  
  A <- pcaplotRK_publication(TSE,intgroup = c("RepsSpecs"), pcX = 1, pcY = 2, title="", ellipse = F, ellipse.prob = 0.5, group_order = COL) +
    scale_fill_manual(values  = COL) +
    scale_color_manual(values = COL) + 
    atheme +
  theme(
         panel.background = element_rect(fill='transparent'), #transparent panel bg
         plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
         ) +
      theme(axis.title.x.bottom =  element_text(color="grey13"), 
        strip.text = element_text(color = "black", face= "bold")) +
      theme(
        legend.title=element_text(size=6), 
        legend.text=element_text(size=6)) +
    theme(
        panel.grid.major = element_line(colour = "grey50"), 
        panel.grid.minor = element_line(colour = "grey50"))

  SampleDist_PCA <- cowplot::plot_grid(A, labels = c(""), ncol = 1)
  
  A <- plot_grid(SampleDist_PCA, ncol = 1)
   plot(A)
  ggsave(A, filename = paste(save_name, gg, "clrPCA.png"), path = pathPlots, device='png', dpi=300, width = 8,
  height = 7)
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
}
## [1] "TSEc.l.r_SSU"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."

## [1] "TSEc.l.r_WF"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."

## [1] "TSEc.l.r_Fish"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."

## [1] "TSEc.l.r_OE"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."

## [1] "TSEc.l.r_GC"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."

## [1] "TSEc.l.r_OEWF"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."

## [1] "TSEc.l.r_GCWF"
## [1] "Adapted from Marini F, Binder H (2019). “pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components.” BMC Bioinformatics, 20(1), 331. doi:10.1186/s12859-019-2879-1, https://bioconductor.org/packages/pcaExplorer/."

6.2 CCA

###################
#CCA with relabund#    
###################
for (Dataset in names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE|ps_Fish", names(pslist)[!grepl("WF", names(pslist))])]){
  require(plyr)
  require(ggrepel) 
  require(cowplot)
  require(DESeq2)
  require(matrixStats)
  tryCatch({
    
  ps <- pslist[[Dataset]]
  
  Datasetname <- sub("ps_", "", Dataset)
  

  
  if (Datasetname == "Fish") {
  COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in% sample_data(ps)$Season]
  GroupingVariable <- c("Season", "LOC")
  } else {
  COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in% sample_data(ps)$Season]
  GroupingVariable <- c("Season", "LOC")
  }

  if (Datasetname %in% c("GC", "Fish")) {
  ps <- prune_samples(sample_names(ps)[!grepl("GCSU21BBEB6", sample_names(ps))], ps) 
  } else {ps <- ps}


  if (Datasetname == "GC") {
  Traits <- c(
    "Salinity","O2","SecchiDepth", "Temperature",
    "Age", "FCF", "HSI", "SSI",
  GroupingVariable)
  print(Traits)
  } else if (Datasetname %in% c("OE", "Fish")) {
  Traits <- c(
    "Salinity","O2","SecchiDepth", "Temperature",
    "Cortison", GroupingVariable)
  print(Traits)
  } else {print("error")}
  
  
  print(paste("Samples removed for low frequency below 10000 seqs in;", g, sep = ""))
  print(which(rowSums(otu_table(ps)) < 5000))
  ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
  ps_Filter_relAbund <- ps_Filter %>%
        transform_sample_counts(function(x) {x/sum(x)*100}) 
  
  sample_data(ps_Filter)<-sample_data(ps_Filter)[ ,(names(sample_data(ps_Filter)) %in% Traits)]
  '%!in%' <- function(x,y)!('%in%'(x,y))
  print("NAs in TraitData")
  print(table(is.na(sample_data(ps_Filter))))

  include <- Traits[Traits != GroupingVariables]

  
  '%!in%' <- function(x,y)!('%in%'(x,y))
  Include <- Traits[Traits %!in% GroupingVariables]

  df_Traits <- data.frame(sample_data(ps_Filter))
  df_Traits <- df_Traits[names(df_Traits) %in% Traits]
  
  df_Include <- df_Traits[names(df_Traits) %in% Include]
  df_Include <- decostand(df_Include, method='standardize')

  ##################################################
  #Determine which variables are useful for the plot#
  ##################################################
  relAbund_table <- as.data.frame(otu_table(ps_Filter_relAbund, taxa_are_rows = FALSE)) 
  
  ##################################################
  #Determine which variables are useful for the plot#
  ##################################################
  abund_table.adonis <- vegan::adonis2(relAbund_table ~ ., data=df_Include)
  print(Datasetname)
  print(abund_table.adonis)
  abund_table.adonis.sig <- abund_table.adonis[!is.na(abund_table.adonis$`Pr(>F)`),]
  abund_table.adonis.sig <- abund_table.adonis.sig[abund_table.adonis.sig$`Pr(>F)` < 0.05,]
  sig.variables <- row.names(abund_table.adonis.sig)
  
  df_Significant <- df_Include[names(df_Include) %in% c(sig.variables)]

  #########################
  #Peform the CCA analysis#
  #########################
  perform_cca <- function(abund_table, meta_table) {
  formula_text <- reformulate(response = "abund_table", 
                            termlabels = names(meta_table))
    tryCatch({
      vegan.cca <- vegan::cca(formula_text, data = meta_table)
      return(vegan.cca)
      }, error = function(e) {
      message("Error performing CCA: ", e$message)
      return(NULL)
      })
  }
  
  result.cca <- perform_cca(relAbund_table, df_Significant)
  scores.cca<- vegan::scores(result.cca, display=c("sp","wa","lc","bp","cn"))
  
  #########################
  #Extract site data first#
  df_sites <- data.frame(scores.cca$sites,df_Traits)

  ############
  #Draw sites#
  p <-  ggplot(data=df_sites,aes(CCA1,CCA2,colour=Season, shape=LOC)) +
        geom_point() +
    scale_color_manual(values = COL) + 
        scale_fill_manual(values = COL) +
    #ggrepel::geom_label_repel(mapping = aes_string(label = "LOC", 
    #        fill = "LOC"), color = "white", show.legend = TRUE)
    ggrepel::geom_text_repel(aes(label = LOC))
    
  multiplier <- vegan:::ordiArrowMul(scores.cca$biplot)
  df_arrows<- scores.cca$biplot*multiplier
  df_arrows=as.data.frame(df_arrows)
  
    p <- p + geom_segment(data=df_arrows,inherit.aes=F, mapping=(aes(x = 0, y = 0, xend = CCA1, yend = CCA2)), arrow = arrow(length = unit(0.2, "cm")),color="grey60",alpha=0.8) +
    geom_text_repel(data=as.data.frame(df_arrows*1.2), inherit.aes=F, mapping=aes(CCA1, CCA2, label =
                    rownames(df_arrows)),color="grey30",alpha=0.99, size = 6)
    
  p <- p + atheme +
   theme(
         panel.background = element_rect(fill='transparent'), #transparent panel bg
         plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
         ) +
      theme(axis.title.x.bottom =  element_text(color="grey13"), 
        strip.text = element_text(color = "black", face= "bold")) +
      theme(
        legend.title=element_text(size=6), 
        legend.text=element_text(size=6)) +
    theme(
        panel.grid.major = element_line(colour = "white"), 
        panel.grid.minor = element_line(colour = "white"))
  
  p <-p+ xlab(paste(sub("CAP1", "", ord_labels(result.cca)[1]))) + 
         ylab(paste(sub("CAP2", "", ord_labels(result.cca)[2]))) 
  
   CCA_plot <- cowplot::plot_grid(p, labels = c(""), ncol = 1)
  
   CCA_plot <- plot_grid(CCA_plot, ncol = 1)
   plot(CCA_plot)
   ggsave(CCA_plot, filename = paste(paste(save_name, "relabund_sigvariable_CCA", Datasetname, sep="_"), ".png", sep=""), path = pathPlots, device='png', dpi=300, width = 8,
   height = 7)
   
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
}
## ERROR : error in evaluating the argument 'table' in selecting a method for function '%in%': could not find function "sample_data" 
## ERROR : error in evaluating the argument 'table' in selecting a method for function '%in%': could not find function "sample_data" 
## ERROR : error in evaluating the argument 'table' in selecting a method for function '%in%': could not find function "sample_data"

6.3 Bray-Dist

We create a average distance dataset here by permutated subsampling

Bray_list <- list()
for (Dataset in c(names(pslist)[grepl("ps_GC|ps_OE|ps_Fish|ps_SSU", names(pslist))], "ps_WF")) {
  
  require(phyloseq)
  require(vegan)
  require(plyr)
  require(dplyr)
  require(ggrepel) 
  require(cowplot)
  require(ggordiplots)
  
  Datasetname <- sub("ps_", "", Dataset)
  #ps <- pslist[[Dataset]]
  
  Samples <- sample_names(pslist[[Dataset]])
  ps <- merge_phyloseq(pslist$ps_OE, pslist$ps_GC, pslist$ps_WF)
  ps <- prune_samples(sample_names(ps) %in% Samples, ps)
  
  
  Bray_list_length <- length(Bray_list)
  
    if (Datasetname %in% c("GC", "OE", "GCWF", "OEWF", "Fish")) {
      print(paste("Samples removed for low frequency below 5000 seqs in;", Dataset, sep = ""))
      print(which(rowSums(otu_table(ps)) < 5000))
      ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
    } else if (Datasetname %in% c("SSU")) {
      print(paste("Samples removed for low frequency below 5000 seqs in;", Dataset, sep = ""))
      print(which(rowSums(otu_table(ps)) < 5000))
      ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
    } else if (Datasetname %in% c("WF")) {
      print(paste("No Waterfiler Samples removed ;", Dataset, sep = ""))
      print(rowSums(otu_table(ps)))
      ps_Filter <- ps
    }
  
  ps_Filter_relAbund <- ps_Filter %>%
        transform_sample_counts(function(x) {x/sum(x)*100}) 
    
  ###############
  #RelAbund Dist#
  ###############
  relAbund_dist_matrix <- 
        as.data.frame(otu_table(ps_Filter_relAbund, taxa_are_rows = FALSE)) %>%
        vegdist(., method="bray")
  
  ##############
  #Average Dist#
  ##############
  smallest_group <- min(rowSums(otu_table(ps_Filter)))
  print("Smallest group for avgdist")
  print(smallest_group)
  Avg_dist_matrix <- avgdist(as.data.frame(otu_table(ps_Filter, 
                        taxa_are_rows = FALSE)), dmethod = "bray", sample = smallest_group)
  
  
  relAbund_dist_dist_tibble <- relAbund_dist_matrix %>%
    as.matrix() %>%
    as_tibble(rownames="sample") %>%
    pivot_longer(-sample) %>%
    filter(name < sample )

  Avg_dist_dist_tibble <- Avg_dist_matrix %>%
    as.matrix() %>%
    as_tibble(rownames="sample") %>%
    pivot_longer(-sample) %>%
    filter(name < sample)
  
  Bray_list[[Bray_list_length+1]] <- relAbund_dist_matrix
  names(Bray_list)[[Bray_list_length+1]] <- paste(Datasetname, "relAbund_dist_matrix", sep="_")

  Bray_list[[Bray_list_length+2]] <- Avg_dist_matrix
  names(Bray_list)[[Bray_list_length+2]] <- paste(Datasetname, "Avg_dist_matrix", sep="_")
  
  Bray_list[[Bray_list_length+3]] <- ps_Filter
  names(Bray_list)[[Bray_list_length+3]] <- paste("ps_Filter", Datasetname, sep="_")
  
}
## Loading required package: phyloseq
## 
## Attaching package: 'phyloseq'
## The following object is masked from 'package:SummarizedExperiment':
## 
##     distance
## The following object is masked from 'package:Biobase':
## 
##     sampleNames
## The following object is masked from 'package:GenomicRanges':
## 
##     distance
## The following object is masked from 'package:IRanges':
## 
##     distance
## Loading required package: vegan
## Warning: package 'vegan' was built under R version 4.3.3
## Loading required package: permute
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 4.3.1
## This is vegan 2.6-6.1
## Loading required package: ggordiplots
## Warning: package 'ggordiplots' was built under R version 4.3.1
## Loading required package: glue
## Warning: package 'glue' was built under R version 4.3.1
## 
## Attaching package: 'glue'
## The following object is masked from 'package:SummarizedExperiment':
## 
##     trim
## The following object is masked from 'package:GenomicRanges':
## 
##     trim
## The following object is masked from 'package:IRanges':
## 
##     trim
## [1] "Samples removed for low frequency below 5000 seqs in;ps_SSU"
## GCAU21TWFL2 GCAU21MLEB1  WFSU21SSFL  WFWI22MGEB 
##         183         184         229         232 
## [1] "Smallest group for avgdist"
## [1] 5998
## [1] "Samples removed for low frequency below 5000 seqs in;ps_Fish"
## GCAU21TWFL2 GCAU21MLEB1 
##         183         184 
## [1] "Smallest group for avgdist"
## [1] 5998
## [1] "Samples removed for low frequency below 5000 seqs in;ps_OE"
## named integer(0)
## [1] "Smallest group for avgdist"
## [1] 9071
## [1] "Samples removed for low frequency below 5000 seqs in;ps_GC"
## GCAU21TWFL2 GCAU21MLEB1 
##          45          46 
## [1] "Smallest group for avgdist"
## [1] 5998
## [1] "Samples removed for low frequency below 5000 seqs in;ps_OEWF"
## WFSU21SSFL WFWI22MGEB 
##        141        144 
## [1] "Smallest group for avgdist"
## [1] 7739
## [1] "Samples removed for low frequency below 5000 seqs in;ps_GCWF"
## GCAU21TWFL2 GCAU21MLEB1  WFSU21SSFL  WFWI22MGEB 
##          45          46          91          94 
## [1] "Smallest group for avgdist"
## [1] 5998
## [1] "No Waterfiler Samples removed ;ps_WF"
## WFSU21MGEB WFSU21BBEB WFSU21SSFL WFSU21MLFL WFAU21TWEB WFWI22MGEB WFWI22MGFL 
##       8983       8409       4572      17204      13785       3576      13431 
## WFWI22BBEB WFWI22BBFL WFWI22MLFL WFWI22SSFL WFSP22MGEB WFSP22MGFL WFSP22BBEB 
##       7739      26491       9265      24818      11908      31196      28110 
## WFSP22BBFL WFSP22SSEB WFSP22SSFL WFSP22MLEB WFSP22MLFL 
##      22786      17448      21514       9152      20503 
## [1] "Smallest group for avgdist"
## [1] 3576
pslist$BrayCurtis <- Bray_list

6.3.1 PCOA

PCOA_list <- list()
  for (Dataset in names(pslist$BrayCurtis)[grepl(
    "OE_Avg_dist_matrix|GC_Avg_dist_matrix|SSU_Avg_dist_matrix|OE_relAbund_dist_matrix|GC_relAbund_dist_matrix|SSU_relAbund_dist_matrix|OEWF_Avg_dist_matrix|GCWF_Avg_dist_matrix", names(pslist$BrayCurtis))]) {
  require(vegan)
  require(plyr)
  require(dplyr)
  require(ggrepel) 
  require(cowplot)
  
  Datasetname      <- sub("_Avg_dist_matrix|_relAbund_dist_matrix", "", Dataset)

  if (grepl("Avg_dist_matrix", Dataset)) {
    BrayKind <- "avgdist" 
    } else if (grepl("relAbund_dist_matrix", Dataset)) {
    BrayKind <- "relabund_dist"
    }
  
  Bray_list_species <- Bray_list[[Dataset]]
  
  PCOA_list_length <- length(PCOA_list)

  COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in%
                                          labels(Bray_list_species)]
  set.seed(1)
  pcoa <- cmdscale(Bray_list_species, eig=TRUE, add=TRUE)
  positions <- pcoa$points
  colnames(positions) <- c("pcoa1", "pcoa2")
  percent_explained <- 100 * pcoa$eig / sum(pcoa$eig)
  
  pretty_pe <- format(round(percent_explained, digits =1), nsmall=1, trim=TRUE)

  labels <- c(glue("PCo Axis 1 ({pretty_pe[1]}%)"),
            glue("PCo Axis 2 ({pretty_pe[2]}%)"))


  if (Datasetname == "SSU") {
    # p <- positions %>%
    #   as_tibble(rownames = "SampleID") %>%
    #   left_join(SAMDF16S, by = c("SampleID" = "SampleID" )) %>%
    #   #left_join(data.frame(sample_data(Bray_list$ps_Filter_OE)), by = c("SampleID" = "SampleID" )) %>%
    #   ggplot(., aes(x=pcoa1, y=pcoa2, color=fSeason, shape=fSpecies)) +
    #   geom_point(size=4) +
    #   labs(x=labels[1], y=labels[2]) + 
    #   #theme_bw() + 
    #   theme(legend.position = "right") +
    #   guides(colour = guide_legend(title.position="top", title.hjust = 0, title = "Season"),
    #      shape = guide_legend(title.position="top", title.hjust = 0, title = "Location")) +
    #   scale_color_manual(values = col.Palette$col.Palette.Season) + 
    #   scale_fill_manual(values = col.Palette$col.Palette.Season) +
    #   ggrepel::geom_text_repel(aes(label = fSpecies)) +
    #   coord_fixed(sqrt(percent_explained[2] / percent_explained[1]))
      
      
  p <- positions %>%
    as_tibble(rownames = "SampleID") %>%
    left_join(SAMDF16S, by = c("SampleID" = "SampleID" )) %>%
    ggplot(aes(x=pcoa1, y=pcoa2,  shape=fSpecies)) +
    geom_point(aes(color=fSeason), size=4) +
    labs(x=labels[1], y=labels[2]) + 
    theme(legend.position = "right") +
    guides(colour = guide_legend(title.position="top", title.hjust = 0, title = "Season"),
         shape = guide_legend(title.position="top", title.hjust = 0, title = "Species")) +
    ggrepel::geom_text_repel(aes(label = fSpecies)) +
    coord_fixed(sqrt(percent_explained[2] / percent_explained[1])) +
    stat_ellipse(aes(fill=fSpecies), geom="polygon", alpha=0.25, colour = "black") +
    scale_color_manual(values = col.Palette$col.Palette.Season) +
    scale_fill_manual(values = col.Palette$col.Palette.Species)
  

    
    
    } else if (Datasetname %in% c("GCWF", "OEWF")) {
  
    p <- positions %>%
    as_tibble(rownames = "SampleID") %>%
    left_join(SAMDF16S, by = c("SampleID" = "SampleID" )) %>%
    ggplot(., aes(x=pcoa1, y=pcoa2)) +
    geom_point(size=4, aes(color=fSeason, shape=fLOC)) +
    labs(x=labels[1], y=labels[2]) + 
    theme(legend.position = "right") +
    #guides(colour = guide_legend(title.position="top", title.hjust = 0, title = "Season"),
    #     shape = guide_legend(title.position="top", title.hjust = 0, title = "Species")) +
    #ggrepel::geom_text_repel(aes(label = SampleID)) +
    coord_fixed(sqrt(percent_explained[2] / percent_explained[1])) +
    stat_ellipse(aes(fill=fSpecies), geom="polygon", alpha=0.25, colour = "black") +
    scale_color_manual(values = col.Palette$col.Palette.Season) + 
    scale_fill_manual(values = col.Palette$col.Palette.Species) +
    scale_shape_manual(values=c(15, 16, 17, 8, 25)) +
      #ggrepel::geom_text_repel(aes(label = fSpecies))+
    coord_fixed(sqrt(percent_explained[2] / percent_explained[1]))
    
    } else {
     p <- positions %>%
      as_tibble(rownames = "SampleID") %>%
      left_join(SAMDF16S, by = c("SampleID" = "SampleID" )) %>%
      #left_join(data.frame(sample_data(Bray_list$ps_Filter_OE)), by = c("SampleID" = "SampleID" )) %>%
      ggplot(., aes(x=pcoa1, y=pcoa2, color=fSeason, shape=fLOC)) +
      geom_point(size=4) +
      labs(x=labels[1], y=labels[2]) + 
      theme_bw()+ theme(legend.position = "right") +
      guides(colour = guide_legend(title.position="top", title.hjust = 0, title = "Season"),
         shape = guide_legend(title.position="top", title.hjust = 0, title = "Location")) +
      scale_color_manual(values = col.Palette$col.Palette.Season) + 
      scale_fill_manual(values = col.Palette$col.Palette.Season) +
      scale_shape_manual(values=c(15, 16, 17, 8, 25)) +
      ggrepel::geom_text_repel(aes(label = fLOC))+
      coord_fixed(sqrt(percent_explained[2] / percent_explained[1]))
     
  }
  
    #p <- p + annotate("text", x = Inf, y = Inf, hjust = 1, vjust = 1, size = 5,
    #       label = paste("Stress =", round(stress, 3)), 
    #       color = c("darkred", "darkred", "darkred"),  fontface = "bold")
               
  
  #p + theme(panel.grid = element_blank())
  
    p <- p + atheme +
      theme(
         panel.background = element_rect(fill='transparent'), #transparent panel bg
         plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
         ) +
      theme(axis.title.x.bottom =  element_text(color="grey13"), 
        strip.text = element_text(color = "black", face= "bold")) +
      theme(
        legend.title=element_text(size=6), 
        legend.text=element_text(size=6)) +
      theme(
        panel.grid.major = element_line(colour = "white"), 
        panel.grid.minor = element_line(colour = "white"))
    p <- p +  theme(legend.position = "bottom", legend.box = "horizontal")
  

   PCOA_plot <- cowplot::plot_grid(p, labels = c(""), ncol = 1)
  
   PCOA_plot <- plot_grid(PCOA_plot, ncol = 1)
   #plot(PCOA_plot)
   ggsave(PCOA_plot, filename = paste(paste(save_name, BrayKind, "PCOA", Datasetname, sep="_"), ".png", sep=""), 
          path = pathPlots, device='png', dpi=300, width = 8,height = 7)
  
   PCOA_list[[PCOA_list_length+1]] <- PCOA_plot
   names(PCOA_list)[[PCOA_list_length+1]] <- paste("PCOA",BrayKind, Datasetname, sep="_")
  }
## Warning: ggrepel: 191 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 195 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## Warning: ggrepel: 69 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## Warning: ggrepel: 70 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## Warning: ggrepel: 34 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## Warning: ggrepel: 32 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will
## replace the existing one.
PCOA_list$PCOA_avgdist_SSU
## Warning: ggrepel: 212 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

6.3.2 NMDS

NMDS_list <- list()
  for (Dataset in names(pslist$BrayCurtis)[grepl(
    "OE_Avg_dist_matrix|GC_Avg_dist_matrix|SSU_Avg_dist_matrix|OE_relAbund_dist_matrix|GC_relAbund_dist_matrix|SSU_relAbund_dist_matrix", names(pslist$BrayCurtis))]) {
  require(vegan)
  require(plyr)
  require(dplyr)
  require(ggrepel) 
  require(cowplot)
  
  Datasetname      <- sub("_Avg_dist_matrix|_relAbund_dist_matrix", "", Dataset)

  if (grepl("Avg_dist_matrix", Dataset)) {
    BrayKind <- "avgdist" 
    } else if (grepl("relAbund_dist_matrix", Dataset)) {
    BrayKind <- "relabund_dist"
    }
  
  Bray_list_species <- Bray_list[[Dataset]]
  
  NMDS_list_length <- length(NMDS_list)

  COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in%
                                          labels(Bray_list_species)]
  set.seed(1)
  nmds <- metaMDS(Bray_list_species)
  stress <- nmds$stress
  
  if (Datasetname == "SSU") {
    p <- scores(nmds) %>%
      as_tibble(rownames = "samples") %>%
      left_join(SAMDF16S, by = c("samples" = "SampleID" )) %>%
      ggplot(aes(x=NMDS1, y=NMDS2, color=fSeason, shape=fLOC)) +
      geom_point(size=4)+ 
      scale_color_manual(values = col.Palette$col.Palette.Season) + 
      scale_fill_manual(values = col.Palette$col.Palette.Season) +
      ggrepel::geom_text_repel(aes(label = fSpecies))
  } else {
    p <- scores(nmds) %>%
      as_tibble(rownames = "samples") %>%
      left_join(SAMDF16S, by = c("samples" = "SampleID" )) %>%
      ggplot(aes(x=NMDS1, y=NMDS2, color=fSeason, shape=fLOC)) +
      geom_point(size=4)+ 
      scale_color_manual(values = col.Palette$col.Palette.Season) + 
      scale_fill_manual(values = col.Palette$col.Palette.Season) +
      scale_shape_manual(values=c(15, 16, 17, 8, 25)) +
      ggrepel::geom_text_repel(aes(label = fLOC))
  }
  
    p <- p + annotate("text", x = Inf, y = Inf, hjust = 1, vjust = 1, size = 5,
           label = paste("Stress =", round(stress, 3)), 
           color = c("darkred", "darkred", "darkred"),  fontface = "bold")
               
    p <- p + atheme +
      theme(
         panel.background = element_rect(fill='transparent'), #transparent panel bg
         plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
         ) +
      theme(axis.title.x.bottom =  element_text(color="grey13"), 
        strip.text = element_text(color = "black", face= "bold")) +
      theme(
        legend.title=element_text(size=6), 
        legend.text=element_text(size=6)) +
      theme(
        panel.grid.major = element_line(colour = "white"), 
        panel.grid.minor = element_line(colour = "white"))
  

   NMDS_plot <- cowplot::plot_grid(p, labels = c(""), ncol = 1)
  
   NMDS_plot <- plot_grid(NMDS_plot, ncol = 1)
   #plot(NMDS_plot)
   ggsave(NMDS_plot, filename = paste(paste(save_name, BrayKind, "NMDS", Datasetname, sep="_"), ".png", sep=""), 
          path = pathPlots, device='png', dpi=300, width = 8,height = 7)
  
   NMDS_list[[NMDS_list_length+1]] <- NMDS_plot
   names(NMDS_list)[[NMDS_list_length+1]] <- paste("NMDS",BrayKind, Datasetname, sep="_")
  }
## Run 0 stress 0.1598705 
## Run 1 stress 0.1589304 
## ... New best solution
## ... Procrustes: rmse 0.01986886  max resid 0.2209935 
## Run 2 stress 0.1603834 
## Run 3 stress 0.1673076 
## Run 4 stress 0.157282 
## ... New best solution
## ... Procrustes: rmse 0.03117025  max resid 0.2195268 
## Run 5 stress 0.1872517 
## Run 6 stress 0.1598648 
## Run 7 stress 0.1587316 
## Run 8 stress 0.1568198 
## ... New best solution
## ... Procrustes: rmse 0.01454307  max resid 0.178522 
## Run 9 stress 0.1632134 
## Run 10 stress 0.1577244 
## Run 11 stress 0.1629727 
## Run 12 stress 0.162905 
## Run 13 stress 0.1625196 
## Run 14 stress 0.1646227 
## Run 15 stress 0.1600792 
## Run 16 stress 0.1625054 
## Run 17 stress 0.1604759 
## Run 18 stress 0.1706289 
## Run 19 stress 0.1629655 
## Run 20 stress 0.1606143 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      3: no. of iterations >= maxit
##     17: stress ratio > sratmax
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## Warning: ggrepel: 156 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Run 0 stress 0.158884 
## Run 1 stress 0.157863 
## ... New best solution
## ... Procrustes: rmse 0.0201235  max resid 0.2197741 
## Run 2 stress 0.1598936 
## Run 3 stress 0.1662939 
## Run 4 stress 0.1573704 
## ... New best solution
## ... Procrustes: rmse 0.02908927  max resid 0.2207995 
## Run 5 stress 0.4174065 
## Run 6 stress 0.1560476 
## ... New best solution
## ... Procrustes: rmse 0.02829525  max resid 0.2234037 
## Run 7 stress 0.1566534 
## Run 8 stress 0.1579184 
## Run 9 stress 0.1622481 
## Run 10 stress 0.1569373 
## Run 11 stress 0.1642362 
## Run 12 stress 0.1629476 
## Run 13 stress 0.1613905 
## Run 14 stress 0.1640273 
## Run 15 stress 0.159106 
## Run 16 stress 0.1614624 
## Run 17 stress 0.1594659 
## Run 18 stress 0.1681732 
## Run 19 stress 0.1615566 
## Run 20 stress 0.1595685 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      2: no. of iterations >= maxit
##     15: stress ratio > sratmax
##      3: scale factor of the gradient < sfgrmin
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## Warning: ggrepel: 162 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Run 0 stress 0.1406636 
## Run 1 stress 0.1642062 
## Run 2 stress 0.1614026 
## Run 3 stress 0.1697188 
## Run 4 stress 0.1544251 
## Run 5 stress 0.1627928 
## Run 6 stress 0.1698772 
## Run 7 stress 0.1398245 
## ... New best solution
## ... Procrustes: rmse 0.01506232  max resid 0.1027648 
## Run 8 stress 0.1541418 
## Run 9 stress 0.1609619 
## Run 10 stress 0.1570996 
## Run 11 stress 0.1423337 
## Run 12 stress 0.1454004 
## Run 13 stress 0.1562743 
## Run 14 stress 0.157439 
## Run 15 stress 0.1666092 
## Run 16 stress 0.1434355 
## Run 17 stress 0.1397957 
## ... New best solution
## ... Procrustes: rmse 0.001201193  max resid 0.01112539 
## Run 18 stress 0.1738992 
## Run 19 stress 0.1426277 
## Run 20 stress 0.1657134 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     13: stress ratio > sratmax
##      7: scale factor of the gradient < sfgrmin
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## Warning: ggrepel: 74 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Run 0 stress 0.1400518 
## Run 1 stress 0.1623834 
## Run 2 stress 0.1605131 
## Run 3 stress 0.1687144 
## Run 4 stress 0.1530577 
## Run 5 stress 0.1626125 
## Run 6 stress 0.1692477 
## Run 7 stress 0.1391221 
## ... New best solution
## ... Procrustes: rmse 0.01506125  max resid 0.1005011 
## Run 8 stress 0.1536073 
## Run 9 stress 0.1598388 
## Run 10 stress 0.1560175 
## Run 11 stress 0.1414518 
## Run 12 stress 0.14463 
## Run 13 stress 0.1556535 
## Run 14 stress 0.1575412 
## Run 15 stress 0.1662978 
## Run 16 stress 0.1426365 
## Run 17 stress 0.1391944 
## ... Procrustes: rmse 0.001685777  max resid 0.01933529 
## Run 18 stress 0.1720325 
## Run 19 stress 0.141849 
## Run 20 stress 0.1606002 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     16: stress ratio > sratmax
##      4: scale factor of the gradient < sfgrmin
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## Warning: ggrepel: 73 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Run 0 stress 0.1636024 
## Run 1 stress 0.1546993 
## ... New best solution
## ... Procrustes: rmse 0.07776629  max resid 0.6089917 
## Run 2 stress 0.1672091 
## Run 3 stress 0.1602318 
## Run 4 stress 0.1589382 
## Run 5 stress 0.156329 
## Run 6 stress 0.1682938 
## Run 7 stress 0.1715344 
## Run 8 stress 0.1740996 
## Run 9 stress 0.171455 
## Run 10 stress 0.1590428 
## Run 11 stress 0.1673151 
## Run 12 stress 0.1619349 
## Run 13 stress 0.1730834 
## Run 14 stress 0.1614785 
## Run 15 stress 0.1594615 
## Run 16 stress 0.154597 
## ... New best solution
## ... Procrustes: rmse 0.01151229  max resid 0.09751585 
## Run 17 stress 0.155883 
## Run 18 stress 0.1698926 
## Run 19 stress 0.1602295 
## Run 20 stress 0.1851819 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      5: no. of iterations >= maxit
##     15: stress ratio > sratmax
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## Warning: ggrepel: 51 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Run 0 stress 0.1638112 
## Run 1 stress 0.1549974 
## ... New best solution
## ... Procrustes: rmse 0.07744348  max resid 0.6090912 
## Run 2 stress 0.1686948 
## Run 3 stress 0.1623907 
## Run 4 stress 0.1586151 
## Run 5 stress 0.1630234 
## Run 6 stress 0.172169 
## Run 7 stress 0.1721439 
## Run 8 stress 0.1677825 
## Run 9 stress 0.1682698 
## Run 10 stress 0.1653278 
## Run 11 stress 0.1675991 
## Run 12 stress 0.161257 
## Run 13 stress 0.1747388 
## Run 14 stress 0.1596532 
## Run 15 stress 0.1672545 
## Run 16 stress 0.1683737 
## Run 17 stress 0.1545751 
## ... New best solution
## ... Procrustes: rmse 0.01840541  max resid 0.1318633 
## Run 18 stress 0.1703611 
## Run 19 stress 0.1622102 
## Run 20 stress 0.1855765 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##      4: no. of iterations >= maxit
##     16: stress ratio > sratmax
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## ggrepel: 51 unlabeled data points (too many overlaps). Consider increasing max.overlaps

6.3.3 PERMANOVA

https://www.nicholas-ollberding.com/post/introduction-to-the-statistical-analysis-of-microbiome-data-in-r/ While PCA is an exploratory data visualization tool, we can test whether the samples cluster beyond that expected by sampling variability using permutational multivariate analysis of variance (PERMANOVA). It does this by partitioning the sums of squares for the within- and between-cluster components using the concept of centroids. Many permutations of the data (i.e. random shuffling) are used to generate the null distribution. The test from ADONIS can be confounded by differences in dispersion (or spread)…so we want to check this as well.

  #########################################################
  #Test null hypotheses for non-rarified and rarified data#
  #########################################################

  # treatment_sample <- sample(treatment_size)
  # adonis2(Bray_list[[paste(Datasetname, "relAbund_dist_matrix", sep="_")]]~ treatment_sample)$`Pr(>F)`[1]
  # #0.269 
  # adonis2(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]]~ treatment_sample)$`Pr(>F)`[1]
  # #0.275
  
  #####################
  #Fish vs Waterfilter#
  #####################
PERMANOVA_list <- list()
  for (Dataset in names(pslist)[grepl("ps_OEWF|ps_GCWF", names(pslist))]){
   require(vegan)
  
    PERMANOVA_list[[Dataset]] <- list()
    Datasetname <- sub("ps_", "", Dataset)
    Bray_list <- pslist$BrayCurtis
    #Other possibility#
    #Generate distance matrix
    #ps_clr <- microbiome::transform(pslist[[Dataset]], "clr")
    #Generate distance matrix
    #clr_dist_matrix <- phyloseq::distance(ps_clr, method = "euclidean") 
    #ADONIS test
    
    print(Datasetname)
    ADONIS <- vegan::adonis2(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Kind)
    print(ADONIS)
    #Dispersion test and plot
    dispr <- vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Kind)
  
    dispr
    plot(dispr, main = "Ordination Centroids and Dispersion Labeled: Aitchison Distance", sub = "")
    boxplot(dispr, main = "", xlab = "")
    vegan::permutest(dispr)
    BetaDisper <- anova(vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Kind))
    #Cant trust Adonis for both
    
   if (BetaDisper$`Pr(>F)`[1] > 0.05) {
     print(Datasetname)
     print("Trust Adonis")
   } else {
     print(Datasetname)
     print("Can not trust Adonis")
   }
    print(BetaDisper )
    print("Next analysis")
    PERMANOVA_list[[Dataset]][[paste("Kind")]]<- ADONIS
  }
## [1] "OEWF"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind)
##                                                                            Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind   1
## Residual                                                                  153
## Total                                                                     154
##                                                                           SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind    6.813
## Residual                                                                    36.228
## Total                                                                       43.041
##                                                                                R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 0.15829
## Residual                                                                  0.84171
## Total                                                                     1.00000
##                                                                                F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 28.773
## Residual                                                                        
## Total                                                                           
##                                                                           Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind  0.001
## Residual                                                                        
## Total                                                                           
##                                                                              
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind ***
## Residual                                                                     
## Total                                                                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "OEWF"
## [1] "Trust Adonis"
## Analysis of Variance Table
## 
## Response: Distances
##            Df Sum Sq   Mean Sq F value Pr(>F)
## Groups      1 0.0022 0.0021677  0.0768 0.7821
## Residuals 153 4.3186 0.0282261               
## [1] "Next analysis"
## [1] "GCWF"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind)
##                                                                            Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind   1
## Residual                                                                  101
## Total                                                                     102
##                                                                           SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind   7.0305
## Residual                                                                   20.5523
## Total                                                                      27.5828
##                                                                                R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 0.25489
## Residual                                                                  0.74511
## Total                                                                     1.00000
##                                                                               F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind 34.55
## Residual                                                                       
## Total                                                                          
##                                                                           Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind  0.001
## Residual                                                                        
## Total                                                                           
##                                                                              
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Kind ***
## Residual                                                                     
## Total                                                                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "GCWF"
## [1] "Trust Adonis"
## Analysis of Variance Table
## 
## Response: Distances
##            Df  Sum Sq  Mean Sq F value Pr(>F)
## Groups      1 0.01827 0.018273  0.7806  0.379
## Residuals 101 2.36418 0.023408               
## [1] "Next analysis"
  ##############################
  #Fish Seasonal vs Waterfilter#
  ##############################
  for (Dataset in names(pslist)[grepl("ps_OEWF|ps_GCWF", names(pslist))]){
   require(vegan)
  
    PERMANOVA_list[[paste(Dataset, Season, sep="_")]] <- list()
    Datasetname <- sub("ps_", "", Dataset)
    
    #Other possibility#
    #Generate distance matrix
    #ps_clr <- microbiome::transform(pslist[[Dataset]], "clr")
    #Generate distance matrix
    #clr_dist_matrix <- phyloseq::distance(ps_clr, method = "euclidean") 
    #ADONIS test
    
    metadata <- sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])
    
    print(Datasetname)
    ADONIS <- vegan::adonis2(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]] ~ gsub(" \\d+$", "", metadata$Reps))
    print(ADONIS)
    #Dispersion test and plot
    dispr <- vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], gsub(" \\d+$", "", metadata$Reps))
  
    dispr
    plot(dispr, main = "Ordination Centroids and Dispersion Labeled: Aitchison Distance", sub = "")
    boxplot(dispr, main = "", xlab = "")
    vegan::permutest(dispr)
    BetaDisper <- anova(vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], gsub(" \\d+$", "", metadata$Reps)))
    #Cant trust Adonis for both
    
   if (BetaDisper$`Pr(>F)`[1] > 0.05) {
     print(Datasetname)
     print("Trust Adonis")
   } else {
     print(Datasetname)
     print("Can not trust Adonis")
   }
    print(BetaDisper )
    print("Next analysis")
    PERMANOVA_list[[paste(Dataset, Season, sep="_")]][[paste("ADONIS")]] <- list()
    PERMANOVA_list[[paste(Dataset, Season, sep="_")]][[paste("ADONIS")]]<- ADONIS
    
    print(Datasetname)

    metadata$RepsSeason <- gsub(" \\d+$", "", metadata$Reps)
    dist_matrix <- Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]]
    PairwiseADONIS <- pairwiseAdonis::pairwise.adonis2(dist_matrix~RepsSeason, data=data.frame(metadata))       
                        
    PERMANOVA_list[[paste(Dataset, Season, sep="_")]][[paste("Pairwise")]] <- list()
    PERMANOVA_list[[paste(Dataset, Season, sep="_")]][[paste("Pairwise")]]<- PairwiseADONIS
     
    PairwiseADONIS_df <- list()
      for (i in names(PairwiseADONIS[-1])) {
      PairwiseADONIS_df[[i]][1] <- PairwiseADONIS[[i]]$`Pr(>F)`[1]
      PairwiseADONIS_df[[i]][2] <- PairwiseADONIS[[i]]$`F`[1]
      #names(PairwiseADONIS_df)[[i]] <- i
      }

  PairwiseADONIS_df <- list()
  for (i in names(PairwiseADONIS)[-1]) {
    PairwiseADONIS_df[[i]] <- data.frame(
      p_value = PairwiseADONIS[[i]]$`Pr(>F)`[1],
      F_statistic = PairwiseADONIS[[i]]$`F`[1]
    )
  }
  PairwiseADONIS_df <- do.call(rbind, PairwiseADONIS_df)
  print(PairwiseADONIS_df)
}
## [1] "OEWF"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ gsub(" \\d+$", "", metadata$Reps))
##                                      Df SumOfSqs      R2      F Pr(>F)    
## gsub(" \\\\d+$", "", metadata$Reps)   4   18.714 0.43478 28.846  0.001 ***
## Residual                            150   24.328 0.56522                  
## Total                               154   43.041 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "OEWF"
## [1] "Can not trust Adonis"
## Analysis of Variance Table
## 
## Response: Distances
##            Df Sum Sq  Mean Sq F value    Pr(>F)    
## Groups      4 1.0837 0.270915  24.721 9.904e-16 ***
## Residuals 150 1.6439 0.010959                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis"
## [1] "OEWF"
##          p_value F_statistic
## SU_vs_AU   0.001   26.214822
## SU_vs_WI   0.001   14.329574
## SU_vs_SP   0.001    6.215647
## SU_vs_WF   0.001   34.369533
## AU_vs_WI   0.001   34.717334
## AU_vs_SP   0.001   39.381833
## AU_vs_WF   0.001   18.085671
## WI_vs_SP   0.001   19.891102
## WI_vs_WF   0.001   47.848857
## SP_vs_WF   0.001   51.100772
## [1] "GCWF"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ gsub(" \\d+$", "", metadata$Reps))
##                                      Df SumOfSqs     R2      F Pr(>F)    
## gsub(" \\\\d+$", "", metadata$Reps)   4   9.7506 0.3535 13.396  0.001 ***
## Residual                             98  17.8322 0.6465                  
## Total                               102  27.5828 1.0000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "GCWF"
## [1] "Can not trust Adonis"
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value    Pr(>F)    
## Groups     4 0.46918 0.117296  6.0819 0.0002054 ***
## Residuals 98 1.89002 0.019286                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis"
## [1] "GCWF"
##          p_value F_statistic
## SU_vs_AU   0.001    6.165993
## SU_vs_WI   0.001    3.731333
## SU_vs_SP   0.001    7.507221
## SU_vs_WF   0.001   20.103508
## AU_vs_WI   0.001    4.375215
## AU_vs_SP   0.001    3.805386
## AU_vs_WF   0.001   29.802728
## WI_vs_SP   0.001    4.289426
## WI_vs_WF   0.001   16.873265
## SP_vs_WF   0.001   36.695406
  df <- PERMANOVA_list$ps_OEWF_SU_AU_WI_SP$Pairwise
    PairwiseADONIS_df <- list()
  for (i in names(df)[-1]) {
    PairwiseADONIS_df[[i]] <- data.frame(
      p_value = df[[i]]$`Pr(>F)`[1],
      F_statistic = df[[i]]$`F`[1]
    )
  }
  PairwiseADONIS_df <- do.call(rbind, PairwiseADONIS_df)
  PairwiseADONIS_df[grepl("WF", rownames(PairwiseADONIS_df)),]
##          p_value F_statistic
## SU_vs_WF   0.001    34.36953
## AU_vs_WF   0.001    18.08567
## WI_vs_WF   0.001    47.84886
## SP_vs_WF   0.001    51.10077
  df <- PERMANOVA_list$ps_GCWF_SU_AU_WI_SP$Pairwise
    PairwiseADONIS_df <- list()
  for (i in names(df)[-1]) {
    PairwiseADONIS_df[[i]] <- data.frame(
      p_value = df[[i]]$`Pr(>F)`[1],
      F_statistic = df[[i]]$`F`[1]
    )
  }
  PairwiseADONIS_df <- do.call(rbind, PairwiseADONIS_df)
  PairwiseADONIS_df[grepl("WF", rownames(PairwiseADONIS_df)),]
##          p_value F_statistic
## SU_vs_WF   0.001    20.10351
## AU_vs_WF   0.001    29.80273
## WI_vs_WF   0.001    16.87326
## SP_vs_WF   0.001    36.69541
  ##################
  #For Fish vs Fish#
  ##################
  #Pairwise Test
  for (Dataset in names(pslist)[grepl("ps_Fish", names(pslist))]){
   require(vegan)
  
    Datasetname <- sub("ps_", "", Dataset)
    
    #Other Approach
    #Generate distance matrix
    # ps_clr <- microbiome::transform(pslist[[Dataset]], "clr")
    # set.seed(123)
    # clr_dist_matrix <- phyloseq::distance(ps_clr, method ="euclidean")
    
    print(Datasetname)
    ADONIS <- vegan::adonis2(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Species)
    print(ADONIS)

    dispr <- vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Species)
    
    dispr
    plot(dispr, main = "Ordination Centroids and Dispersion Labeled: Aitchison Distance", sub = "")
    boxplot(dispr, main = "", xlab = "")
    vegan::permutest(dispr)
  BetaDisper <- anova(vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Species))
    #Cant trust Adonis for both

   if (BetaDisper$`Pr(>F)`[1] > 0.05) {
     print(Datasetname)
     print("Trust Adonis")
   } else {
     print(Datasetname)
     print("Can not trust Adonis")
   }
    print(BetaDisper )
    print("Next analysis")
    PERMANOVA_list[[Dataset]][[paste("Species")]]<- ADONIS
  
}
## [1] "Fish"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species)
##                                                                               Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species   1
## Residual                                                                     222
## Total                                                                        223
##                                                                              SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species    1.353
## Residual                                                                       49.747
## Total                                                                          51.100
##                                                                                   R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species 0.02647
## Residual                                                                     0.97353
## Total                                                                        1.00000
##                                                                                   F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species 6.0362
## Residual                                                                           
## Total                                                                              
##                                                                              Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species  0.001
## Residual                                                                           
## Total                                                                              
##                                                                                 
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Species ***
## Residual                                                                        
## Total                                                                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "Fish"
## [1] "Trust Adonis"
## Analysis of Variance Table
## 
## Response: Distances
##            Df Sum Sq Mean Sq F value  Pr(>F)  
## Groups      1 0.1038 0.10381  3.6425 0.05761 .
## Residuals 222 6.3270 0.02850                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis"
  ###################
  #For Fish Seasonal#
  ###################
  #Pairwise Test
  for (Dataset in c(names(pslist)[!grepl("WF", names(pslist))][grepl("ps_OE|ps_GC", names(pslist)[!grepl("WF", names(pslist))])], "ps_WF")){
   require(vegan)
    PERMANOVA_list[[Dataset]][[paste("Season")]] <- list()
    Datasetname <- sub("ps_", "", Dataset)
    
    #Other Approach
    #Generate distance matrix
    # ps_clr <- microbiome::transform(pslist[[Dataset]], "clr")
    # set.seed(123)
    # clr_dist_matrix <- phyloseq::distance(ps_clr, method ="euclidean")
    
    print(Datasetname)
    ADONIS <- vegan::adonis2(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Season)
    print(ADONIS)

    dispr <- vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Season)
    
    dispr
    plot(dispr, main = "Ordination Centroids and Dispersion Labeled: Aitchison Distance", sub = "")
    boxplot(dispr, main = "", xlab = "")
    vegan::permutest(dispr)
  BetaDisper <- anova(vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]], sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$Season))
    #Cant trust Adonis for both

   if (BetaDisper$`Pr(>F)`[1] > 0.05) {
     print(Datasetname)
     print("Trust Adonis")
   } else {
     print(Datasetname)
     print("Can not trust Adonis")
   }
    print(BetaDisper )
    print("Next analysis")
    
    PERMANOVA_list[[Dataset]][[paste("Season")]][[paste("Adonis")]] <- list()
    PERMANOVA_list[[Dataset]][[paste("Season")]][[paste("Adonis")]] <- ADONIS
  
    metadata <- sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])
    dist_matrix <- Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]]
    print(Datasetname)
    PairwiseADONIS <- pairwiseAdonis::pairwise.adonis2(dist_matrix~Season,data=data.frame(metadata))       
    PairwiseADONIS                               
    PERMANOVA_list[[Dataset]][[paste("Season")]][[paste("Pairwise")]]<-PairwiseADONIS
     
  }
## [1] "OE"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season)
##                                                                              Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season   3
## Residual                                                                    134
## Total                                                                       137
##                                                                             SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season   11.909
## Residual                                                                      20.664
## Total                                                                         32.573
##                                                                                  R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 0.36562
## Residual                                                                    0.63438
## Total                                                                       1.00000
##                                                                                  F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 25.743
## Residual                                                                          
## Total                                                                             
##                                                                             Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season  0.001
## Residual                                                                          
## Total                                                                             
##                                                                                
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season ***
## Residual                                                                       
## Total                                                                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "OE"
## [1] "Can not trust Adonis"
## Analysis of Variance Table
## 
## Response: Distances
##            Df  Sum Sq Mean Sq F value    Pr(>F)    
## Groups      3 0.98489 0.32830  29.441 1.097e-14 ***
## Residuals 134 1.49422 0.01115                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis"
## [1] "OE"
## [1] "GC"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season)
##                                                                             Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season  3
## Residual                                                                    82
## Total                                                                       85
##                                                                             SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season   2.7221
## Residual                                                                     14.2016
## Total                                                                        16.9237
##                                                                                  R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 0.16084
## Residual                                                                    0.83916
## Total                                                                       1.00000
##                                                                                  F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 5.2391
## Residual                                                                          
## Total                                                                             
##                                                                             Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season  0.001
## Residual                                                                          
## Total                                                                             
##                                                                                
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season ***
## Residual                                                                       
## Total                                                                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "GC"
## [1] "Can not trust Adonis"
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value    Pr(>F)    
## Groups     3 0.38924 0.129746  6.1046 0.0008426 ***
## Residuals 82 1.74281 0.021254                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis"
## [1] "GC"
## [1] "WF"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season)
##                                                                             Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season  3
## Residual                                                                    15
## Total                                                                       18
##                                                                             SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season   1.4831
## Residual                                                                      2.8284
## Total                                                                         4.3115
##                                                                                  R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 0.34398
## Residual                                                                    0.65602
## Total                                                                       1.00000
##                                                                                  F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season 2.6218
## Residual                                                                          
## Total                                                                             
##                                                                             Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season  0.001
## Residual                                                                          
## Total                                                                             
##                                                                                
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$Season ***
## Residual                                                                       
## Total                                                                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "WF"
## [1] "Can not trust Adonis"
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value  Pr(>F)  
## Groups     3 0.16145 0.053816  4.2135 0.02386 *
## Residuals 15 0.19159 0.012772                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis"
## [1] "WF"
  ##################
  #For Fish Spatial#
  ##################
  for (Dataset in names(pslist)[!grepl("WF", names(pslist))][grepl("ps_OE|ps_GC", names(pslist)[!grepl("WF", names(pslist))])]){
   require(vegan)
    PERMANOVA_list[[Dataset]][[paste("Replicates")]] <- list()
    Datasetname <- sub("ps_", "", Dataset)
    
    #Other Approach
    #Generate distance matrix
    # ps_clr <- microbiome::transform(pslist[[Dataset]], "clr")
    # set.seed(123)
    # clr_dist_matrix <- phyloseq::distance(ps_clr, method ="euclidean")
    
    print(Datasetname)
    ADONIS <- vegan::adonis2(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]] ~
              sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$fReplicates)
    print(ADONIS)

    dispr <- vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]],
                    sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$fReplicates)
    
    dispr
    plot(dispr, main = "Ordination Centroids and Dispersion Labeled: Aitchison Distance", sub = "")
    boxplot(dispr, main = "", xlab = "")
    vegan::permutest(dispr)
  BetaDisper <- anova(vegan::betadisper(Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]],
                sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])$fReplicates))
    #Cant trust Adonis for both

   if (BetaDisper$`Pr(>F)`[1] > 0.05) {
     print(Datasetname)
     print("Trust Adonis")
   } else {
     print(Datasetname)
     print("Can not trust Adonis")
   }
    print(BetaDisper )
    print("Next analysis")
    
    PERMANOVA_list[[Dataset]][[paste("Replicates")]][[paste("Adonis")]] <- list()
    PERMANOVA_list[[Dataset]][[paste("Replicates")]][[paste("Adonis")]] <- ADONIS
  
    metadata <- sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep="_")]])
    dist_matrix <- Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep="_")]]
    
    PairwiseADONIS <- pairwiseAdonis::pairwise.adonis2(dist_matrix~fReplicates,data=data.frame(metadata))       
    PERMANOVA_list[[Dataset]][[paste("Replicates")]][[paste("Pairwise")]] <- list()
    PERMANOVA_list[[Dataset]][[paste("Replicates")]][[paste("Pairwise")]]<- PairwiseADONIS
    
  PairwiseADONIS_df <- list()
  for (i in names(PairwiseADONIS[-1])) {
    PairwiseADONIS_df[[i]] <- PairwiseADONIS[[i]]$`Pr(>F)`[1]
    #names(PairwiseADONIS_df)[[i]] <- i
  }

  PairwiseADONIS_df <- PairwiseADONIS_df %>%
    do.call(rbind, .) %>%      
    as.data.frame() %>%   
    tibble::rownames_to_column(var = "Comparison") %>%  
    dplyr::rename(p=2) %>%                       
    filter(p < 0.05)
  print(PairwiseADONIS_df)

}
## [1] "OE"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates)
##                                                                                   Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates  19
## Residual                                                                         118
## Total                                                                            137
##                                                                                  SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates   19.874
## Residual                                                                           12.699
## Total                                                                              32.573
##                                                                                       R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 0.61015
## Residual                                                                         0.38985
## Total                                                                            1.00000
##                                                                                       F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 9.7198
## Residual                                                                               
## Total                                                                                  
##                                                                                  Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates  0.001
## Residual                                                                               
## Total                                                                                  
##                                                                                     
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates ***
## Residual                                                                            
## Total                                                                               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "OE"
## [1] "Can not trust Adonis"
## Analysis of Variance Table
## 
## Response: Distances
##            Df  Sum Sq  Mean Sq F value    Pr(>F)    
## Groups     19 0.76993 0.040523   4.652 7.784e-08 ***
## Residuals 118 1.02789 0.008711                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Next analysis"
##               Comparison     p
## 1   OESU21MG_vs_OESU21SS 0.001
## 2   OESU21MG_vs_OESU21TW 0.002
## 3   OESU21MG_vs_OESU21ML 0.003
## 4   OESU21MG_vs_OEAU21MG 0.001
## 5   OESU21MG_vs_OEAU21BB 0.001
## 6   OESU21MG_vs_OEAU21SS 0.003
## 7   OESU21MG_vs_OEAU21TW 0.001
## 8   OESU21MG_vs_OEAU21ML 0.001
## 9   OESU21MG_vs_OEWI22BB 0.003
## 10  OESU21MG_vs_OEWI22MG 0.001
## 11  OESU21MG_vs_OEWI22SS 0.001
## 12  OESU21MG_vs_OEWI22TW 0.002
## 13  OESU21MG_vs_OEWI22ML 0.003
## 14  OESU21MG_vs_OESP22MG 0.002
## 15  OESU21MG_vs_OESP22BB 0.002
## 16  OESU21MG_vs_OESP22SS 0.001
## 17  OESU21MG_vs_OESP22TW 0.003
## 18  OESU21MG_vs_OESP22ML 0.001
## 19  OESU21BB_vs_OESU21SS 0.017
## 20  OESU21BB_vs_OESU21TW 0.002
## 21  OESU21BB_vs_OESU21ML 0.001
## 22  OESU21BB_vs_OEAU21MG 0.001
## 23  OESU21BB_vs_OEAU21BB 0.001
## 24  OESU21BB_vs_OEAU21SS 0.001
## 25  OESU21BB_vs_OEAU21TW 0.001
## 26  OESU21BB_vs_OEAU21ML 0.003
## 27  OESU21BB_vs_OEWI22BB 0.004
## 28  OESU21BB_vs_OEWI22MG 0.001
## 29  OESU21BB_vs_OEWI22SS 0.001
## 30  OESU21BB_vs_OEWI22TW 0.002
## 31  OESU21BB_vs_OEWI22ML 0.002
## 32  OESU21BB_vs_OESP22MG 0.002
## 33  OESU21BB_vs_OESP22BB 0.038
## 34  OESU21BB_vs_OESP22SS 0.006
## 35  OESU21BB_vs_OESP22TW 0.001
## 36  OESU21BB_vs_OESP22ML 0.011
## 37  OESU21SS_vs_OESU21TW 0.001
## 38  OESU21SS_vs_OESU21ML 0.004
## 39  OESU21SS_vs_OEAU21MG 0.001
## 40  OESU21SS_vs_OEAU21BB 0.001
## 41  OESU21SS_vs_OEAU21SS 0.001
## 42  OESU21SS_vs_OEAU21TW 0.002
## 43  OESU21SS_vs_OEAU21ML 0.002
## 44  OESU21SS_vs_OEWI22BB 0.002
## 45  OESU21SS_vs_OEWI22MG 0.001
## 46  OESU21SS_vs_OEWI22SS 0.002
## 47  OESU21SS_vs_OEWI22TW 0.001
## 48  OESU21SS_vs_OEWI22ML 0.001
## 49  OESU21SS_vs_OESP22MG 0.001
## 50  OESU21SS_vs_OESP22BB 0.001
## 51  OESU21SS_vs_OESP22SS 0.001
## 52  OESU21SS_vs_OESP22TW 0.001
## 53  OESU21SS_vs_OESP22ML 0.008
## 54  OESU21TW_vs_OESU21ML 0.029
## 55  OESU21TW_vs_OEAU21MG 0.001
## 56  OESU21TW_vs_OEAU21BB 0.002
## 57  OESU21TW_vs_OEAU21SS 0.001
## 58  OESU21TW_vs_OEAU21TW 0.002
## 59  OESU21TW_vs_OEAU21ML 0.002
## 60  OESU21TW_vs_OEWI22BB 0.001
## 61  OESU21TW_vs_OEWI22MG 0.001
## 62  OESU21TW_vs_OEWI22SS 0.001
## 63  OESU21TW_vs_OEWI22TW 0.001
## 64  OESU21TW_vs_OEWI22ML 0.001
## 65  OESU21TW_vs_OESP22MG 0.002
## 66  OESU21TW_vs_OESP22BB 0.001
## 67  OESU21TW_vs_OESP22SS 0.001
## 68  OESU21TW_vs_OESP22TW 0.003
## 69  OESU21TW_vs_OESP22ML 0.001
## 70  OESU21ML_vs_OEAU21MG 0.002
## 71  OESU21ML_vs_OEAU21BB 0.002
## 72  OESU21ML_vs_OEAU21SS 0.003
## 73  OESU21ML_vs_OEAU21TW 0.003
## 74  OESU21ML_vs_OEAU21ML 0.002
## 75  OESU21ML_vs_OEWI22BB 0.002
## 76  OESU21ML_vs_OEWI22MG 0.001
## 77  OESU21ML_vs_OEWI22SS 0.001
## 78  OESU21ML_vs_OEWI22TW 0.004
## 79  OESU21ML_vs_OEWI22ML 0.002
## 80  OESU21ML_vs_OESP22MG 0.003
## 81  OESU21ML_vs_OESP22BB 0.002
## 82  OESU21ML_vs_OESP22SS 0.004
## 83  OESU21ML_vs_OESP22TW 0.002
## 84  OESU21ML_vs_OESP22ML 0.004
## 85  OEAU21MG_vs_OEAU21BB 0.001
## 86  OEAU21MG_vs_OEAU21SS 0.001
## 87  OEAU21MG_vs_OEAU21TW 0.002
## 88  OEAU21MG_vs_OEAU21ML 0.002
## 89  OEAU21MG_vs_OEWI22BB 0.001
## 90  OEAU21MG_vs_OEWI22MG 0.002
## 91  OEAU21MG_vs_OEWI22SS 0.001
## 92  OEAU21MG_vs_OEWI22TW 0.004
## 93  OEAU21MG_vs_OEWI22ML 0.002
## 94  OEAU21MG_vs_OESP22MG 0.001
## 95  OEAU21MG_vs_OESP22BB 0.001
## 96  OEAU21MG_vs_OESP22SS 0.006
## 97  OEAU21MG_vs_OESP22TW 0.001
## 98  OEAU21MG_vs_OESP22ML 0.004
## 99  OEAU21BB_vs_OEAU21SS 0.004
## 100 OEAU21BB_vs_OEAU21TW 0.001
## 101 OEAU21BB_vs_OEAU21ML 0.001
## 102 OEAU21BB_vs_OEWI22BB 0.001
## 103 OEAU21BB_vs_OEWI22MG 0.001
## 104 OEAU21BB_vs_OEWI22SS 0.001
## 105 OEAU21BB_vs_OEWI22TW 0.001
## 106 OEAU21BB_vs_OEWI22ML 0.001
## 107 OEAU21BB_vs_OESP22MG 0.002
## 108 OEAU21BB_vs_OESP22BB 0.003
## 109 OEAU21BB_vs_OESP22SS 0.001
## 110 OEAU21BB_vs_OESP22TW 0.002
## 111 OEAU21BB_vs_OESP22ML 0.001
## 112 OEAU21SS_vs_OEAU21TW 0.007
## 113 OEAU21SS_vs_OEAU21ML 0.001
## 114 OEAU21SS_vs_OEWI22BB 0.001
## 115 OEAU21SS_vs_OEWI22MG 0.003
## 116 OEAU21SS_vs_OEWI22SS 0.002
## 117 OEAU21SS_vs_OEWI22TW 0.001
## 118 OEAU21SS_vs_OEWI22ML 0.001
## 119 OEAU21SS_vs_OESP22MG 0.001
## 120 OEAU21SS_vs_OESP22BB 0.001
## 121 OEAU21SS_vs_OESP22SS 0.001
## 122 OEAU21SS_vs_OESP22TW 0.001
## 123 OEAU21SS_vs_OESP22ML 0.001
## 124 OEAU21TW_vs_OEAU21ML 0.008
## 125 OEAU21TW_vs_OEWI22BB 0.001
## 126 OEAU21TW_vs_OEWI22MG 0.001
## 127 OEAU21TW_vs_OEWI22SS 0.001
## 128 OEAU21TW_vs_OEWI22TW 0.001
## 129 OEAU21TW_vs_OEWI22ML 0.003
## 130 OEAU21TW_vs_OESP22MG 0.002
## 131 OEAU21TW_vs_OESP22BB 0.002
## 132 OEAU21TW_vs_OESP22SS 0.001
## 133 OEAU21TW_vs_OESP22TW 0.001
## 134 OEAU21TW_vs_OESP22ML 0.001
## 135 OEAU21ML_vs_OEWI22BB 0.002
## 136 OEAU21ML_vs_OEWI22MG 0.002
## 137 OEAU21ML_vs_OEWI22SS 0.001
## 138 OEAU21ML_vs_OEWI22TW 0.003
## 139 OEAU21ML_vs_OEWI22ML 0.002
## 140 OEAU21ML_vs_OESP22MG 0.001
## 141 OEAU21ML_vs_OESP22BB 0.001
## 142 OEAU21ML_vs_OESP22SS 0.001
## 143 OEAU21ML_vs_OESP22TW 0.001
## 144 OEAU21ML_vs_OESP22ML 0.001
## 145 OEWI22BB_vs_OEWI22TW 0.012
## 146 OEWI22BB_vs_OESP22MG 0.001
## 147 OEWI22BB_vs_OESP22BB 0.001
## 148 OEWI22BB_vs_OESP22SS 0.001
## 149 OEWI22BB_vs_OESP22TW 0.001
## 150 OEWI22BB_vs_OESP22ML 0.001
## 151 OEWI22MG_vs_OEWI22TW 0.027
## 152 OEWI22MG_vs_OESP22MG 0.001
## 153 OEWI22MG_vs_OESP22BB 0.001
## 154 OEWI22MG_vs_OESP22SS 0.001
## 155 OEWI22MG_vs_OESP22TW 0.002
## 156 OEWI22MG_vs_OESP22ML 0.002
## 157 OEWI22SS_vs_OEWI22TW 0.045
## 158 OEWI22SS_vs_OESP22MG 0.001
## 159 OEWI22SS_vs_OESP22BB 0.001
## 160 OEWI22SS_vs_OESP22SS 0.001
## 161 OEWI22SS_vs_OESP22TW 0.001
## 162 OEWI22SS_vs_OESP22ML 0.001
## 163 OEWI22TW_vs_OEWI22ML 0.019
## 164 OEWI22TW_vs_OESP22MG 0.001
## 165 OEWI22TW_vs_OESP22BB 0.001
## 166 OEWI22TW_vs_OESP22SS 0.001
## 167 OEWI22TW_vs_OESP22TW 0.003
## 168 OEWI22TW_vs_OESP22ML 0.001
## 169 OEWI22ML_vs_OESP22MG 0.001
## 170 OEWI22ML_vs_OESP22BB 0.001
## 171 OEWI22ML_vs_OESP22SS 0.003
## 172 OEWI22ML_vs_OESP22TW 0.001
## 173 OEWI22ML_vs_OESP22ML 0.001
## 174 OESP22MG_vs_OESP22BB 0.003
## 175 OESP22MG_vs_OESP22SS 0.001
## 176 OESP22MG_vs_OESP22TW 0.002
## 177 OESP22MG_vs_OESP22ML 0.001
## 178 OESP22BB_vs_OESP22SS 0.001
## 179 OESP22BB_vs_OESP22TW 0.001
## 180 OESP22BB_vs_OESP22ML 0.001
## 181 OESP22SS_vs_OESP22TW 0.002
## 182 OESP22SS_vs_OESP22ML 0.039
## 183 OESP22TW_vs_OESP22ML 0.001
## [1] "GC"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = Bray_list[[paste(Datasetname, "Avg_dist_matrix", sep = "_")]] ~ sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates)
##                                                                                  Df
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 14
## Residual                                                                         71
## Total                                                                            85
##                                                                                  SumOfSqs
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates   6.2966
## Residual                                                                          10.6271
## Total                                                                             16.9237
##                                                                                       R2
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 0.37206
## Residual                                                                         0.62794
## Total                                                                            1.00000
##                                                                                       F
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates 3.0048
## Residual                                                                               
## Total                                                                                  
##                                                                                  Pr(>F)
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates  0.001
## Residual                                                                               
## Total                                                                                  
##                                                                                     
## sample_data(Bray_list[[paste("ps_Filter", Datasetname, sep = "_")]])$fReplicates ***
## Residual                                                                            
## Total                                                                               
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## [1] "GC"
## [1] "Trust Adonis"
## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Groups    14 0.50195 0.035854  1.5555 0.1141
## Residuals 71 1.63652 0.023050               
## [1] "Next analysis"
##              Comparison     p
## 1  GCSU21BB_vs_GCSU21TW 0.008
## 2  GCSU21BB_vs_GCSU21ML 0.001
## 3  GCSU21BB_vs_GCAU21BB 0.022
## 4  GCSU21BB_vs_GCAU21SS 0.005
## 5  GCSU21BB_vs_GCAU21ML 0.037
## 6  GCSU21BB_vs_GCWI22BB 0.035
## 7  GCSU21BB_vs_GCWI22SS 0.005
## 8  GCSU21BB_vs_GCWI22ML 0.016
## 9  GCSU21BB_vs_GCSP22BB 0.032
## 10 GCSU21BB_vs_GCSP22TW 0.044
## 11 GCSU21BB_vs_GCSP22ML 0.007
## 12 GCSU21SS_vs_GCSU21TW 0.036
## 13 GCSU21SS_vs_GCSU21ML 0.005
## 14 GCSU21SS_vs_GCAU21BB 0.012
## 15 GCSU21SS_vs_GCAU21SS 0.001
## 16 GCSU21SS_vs_GCAU21ML 0.014
## 17 GCSU21SS_vs_GCWI22BB 0.040
## 18 GCSU21SS_vs_GCWI22SS 0.013
## 19 GCSU21SS_vs_GCSP22BB 0.002
## 20 GCSU21SS_vs_GCSP22ML 0.001
## 21 GCSU21TW_vs_GCSU21ML 0.030
## 22 GCSU21TW_vs_GCAU21BB 0.004
## 23 GCSU21TW_vs_GCAU21SS 0.002
## 24 GCSU21TW_vs_GCAU21ML 0.004
## 25 GCSU21TW_vs_GCWI22BB 0.006
## 26 GCSU21TW_vs_GCWI22SS 0.001
## 27 GCSU21TW_vs_GCWI22ML 0.012
## 28 GCSU21TW_vs_GCSP22BB 0.003
## 29 GCSU21TW_vs_GCSP22SS 0.005
## 30 GCSU21TW_vs_GCSP22TW 0.006
## 31 GCSU21TW_vs_GCSP22ML 0.001
## 32 GCSU21ML_vs_GCAU21BB 0.001
## 33 GCSU21ML_vs_GCAU21SS 0.001
## 34 GCSU21ML_vs_GCAU21ML 0.001
## 35 GCSU21ML_vs_GCWI22SS 0.003
## 36 GCSU21ML_vs_GCWI22ML 0.014
## 37 GCSU21ML_vs_GCSP22BB 0.001
## 38 GCSU21ML_vs_GCSP22SS 0.001
## 39 GCSU21ML_vs_GCSP22TW 0.001
## 40 GCSU21ML_vs_GCSP22ML 0.001
## 41 GCAU21BB_vs_GCAU21SS 0.006
## 42 GCAU21BB_vs_GCWI22BB 0.037
## 43 GCAU21BB_vs_GCWI22SS 0.019
## 44 GCAU21BB_vs_GCSP22BB 0.036
## 45 GCAU21BB_vs_GCSP22SS 0.024
## 46 GCAU21BB_vs_GCSP22ML 0.001
## 47 GCAU21SS_vs_GCWI22BB 0.021
## 48 GCAU21SS_vs_GCWI22SS 0.002
## 49 GCAU21SS_vs_GCWI22ML 0.005
## 50 GCAU21SS_vs_GCSP22BB 0.006
## 51 GCAU21SS_vs_GCSP22SS 0.012
## 52 GCAU21SS_vs_GCSP22TW 0.033
## 53 GCAU21SS_vs_GCSP22ML 0.002
## 54 GCAU21TW_vs_GCSP22BB 0.030
## 55 GCAU21TW_vs_GCSP22ML 0.032
## 56 GCAU21ML_vs_GCWI22BB 0.041
## 57 GCAU21ML_vs_GCWI22SS 0.044
## 58 GCAU21ML_vs_GCSP22BB 0.012
## 59 GCWI22BB_vs_GCSP22BB 0.022
## 60 GCWI22BB_vs_GCSP22SS 0.032
## 61 GCWI22SS_vs_GCWI22ML 0.034
## 62 GCWI22SS_vs_GCSP22BB 0.001
## 63 GCWI22SS_vs_GCSP22SS 0.038
## 64 GCWI22SS_vs_GCSP22TW 0.028
## 65 GCWI22ML_vs_GCSP22BB 0.009
## 66 GCWI22ML_vs_GCSP22SS 0.030
## 67 GCWI22ML_vs_GCSP22TW 0.030
## 68 GCWI22ML_vs_GCSP22ML 0.011
## 69 GCSP22BB_vs_GCSP22TW 0.034
## 70 GCSP22BB_vs_GCSP22ML 0.005
pslist$PERMANOVA_list <- PERMANOVA_list

https://astrobiomike.github.io/amplicon/dada2_workflow_ex

Checking by all sample types, we get a significant result (0.002) from the betadisper test. This tells us that there is a difference between group dispersions, which means that we can’t trust the results of an adonis (permutational anova) test on this, because the assumption of homogenous within-group disperions is not met. This isn’t all that surprising considering how different the water and biofilm samples are from the rocks.

An important note: An NMDS is just a visualization technique, and is not a statistical assessment of sample separation or correlation. For this you should run a statistical test, like an ANOSIM test for categorical variables, or a Mantel test for continuous variables. Other ordination techniques like PCA, CA, CCA, RDA, etc, may be more useful to you if you want your axes to be meaningful, or if you want to talk about variation partitioning.

6.3.4 dbRDA / CAP

distance based RDA automatically removing Oxygen values for multicollinearity with almost all other abiotic measurements, further there is an automatic exclusion of other multicollinear factors, inspect the output datasets instead of trusting the RDA plot from scratch

RDA_list <- list()
for (Dataset in c(names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE", names(pslist)[!grepl("WF", names(pslist))])], "ps_WF")) {
  
  require(phyloseq)
  require(vegan)
  require(plyr)
  require(dplyr)
  require(ggrepel) 
  require(cowplot)
  require(ggordiplots)
  
  Datasetname <- sub("ps_", "", Dataset)
  ps <- pslist[[Dataset]]
  RDA_list_length <- length(RDA_list)
  
  GroupingVariables <- c("fSeason", "fLOC")
  
  if (Datasetname %in% c("GC", "OE")) {
    
    Traits <- c(
      "Salinity", "Temperature",
      "Age", "FCF", "HSI", "SSI", "GSI", "FillLevel", 
      "NH4",  "NO2", "NO3", "O2", "PO4", "pH", "SPM", "TOC",
      GroupingVariables)
      print(Traits)
    } else if (Datasetname %in% c("WF")) {
    Traits <- c(
      "Salinity", "Temperature", "NH4",  "NO2", "NO3", "O2", "PO4", "pH", "SPM", "TOC", GroupingVariables)
    print(Traits)
  
    } else {print("error")}
  
    if (Datasetname %in% c("GC", "OE")) {
      print(paste("Samples removed for low frequency below 5000 seqs in;", Dataset, sep = ""))
      print(which(rowSums(otu_table(ps)) < 5000))
      ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
    } else if (Datasetname %in% c("WF")) {
      print(paste("No Waterfiler Samples removed ;", Dataset, sep = ""))
      ps_Filter <- ps
    }
    
  '%!in%' <- function(x,y)!('%in%'(x,y))
  Include <- Traits[Traits %!in% GroupingVariables]

  df_Traits <- data.frame(sample_data(ps_Filter))

  df_Traits <- df_Traits[names(df_Traits) %in% Traits]

  df_Include <- df_Traits[names(df_Traits) %in% Include]
  
  print("NAs in TraitData")
  print(table(is.na(df_Include)))
  
  print((df_Include[!complete.cases(df_Include), ]))
  
  df_Include  <- df_Include[complete.cases(df_Include), ]
  
  df_Traits  <- df_Traits[complete.cases(df_Traits), ]
  dim(df_Traits)
  
  df_Include <- decostand(df_Include, method='standardize')
  dim(df_Include)
  
  ps_Filter <- prune_samples(sample_names(ps_Filter)[sample_names(ps_Filter) %in%
               rownames(df_Include[complete.cases(df_Include), ])], ps_Filter)

  ##################################################
  #Determine which variables are useful for the plot#
  ##################################################
  abund_table <- as.data.frame(otu_table(ps_Filter, taxa_are_rows = FALSE)) 
  print("Smallest Group for avgdist")
  smallest_group <- min(rowSums(otu_table(ps_Filter)))
  print(smallest_group)
  avgdist_abund_table  <- avgdist(abund_table, dmethod = "bray", sample = smallest_group)
  dpRDA_results <- capscale(avgdist_abund_table ~ ., data=df_Include, distance='bray')
  
  abund_table.adonis <- vegan::adonis2(avgdist_abund_table ~ ., data=df_Include)
  print(paste("PERMANOVA", Datasetname))
  print(abund_table.adonis)
  abund_table.adonis.sig <- abund_table.adonis[!is.na(abund_table.adonis$`Pr(>F)`),]
  abund_table.adonis.sig <- abund_table.adonis.sig[abund_table.adonis.sig$`Pr(>F)` < 0.05,]
  PERMANOVAsig.variables <- row.names(abund_table.adonis.sig)
  
  ###########################################################################
  #Compute Distance-based Redundancy Analysis based on average bray distance#
  ###########################################################################
  dbRDA_results <- capscale(avgdist_abund_table ~ ., data=df_Include, distance='bray')

  dbrda_plot = plot(dpRDA_results, xlab=ord_labels(dpRDA_results)[1], ylab=ord_labels(dpRDA_results)[2])
  anova(dbRDA_results) # permutation test of statistical significance

  ####################
  #Variable selection#
  ####################
    
    ######################
    #Check for redundancy#
    ######################
    print("VIF")
    print(sort(vif.cca(dpRDA_results))) # calculate VIF, scores >10 are redundant 
    vars.envfit <- names(which(vif.cca(dpRDA_results) <= 10))

    ##################
    # model selection#
    ##################
    rda1 <- capscale(avgdist_abund_table ~ ., data=df_Include) # set up full and null models for 'ordistep' full model
    rda0 <- capscale(avgdist_abund_table ~ 1, data=df_Include) # intercept-only (null) model
    step.env <- ordistep(rda0, scope=formula(rda1), direction='both')  
    #### # perform forward and backward selection of explanatory variables
    print(paste("Ordisep", Datasetname))
    print(step.env$anova)
    vars.ordistep <- gsub('^. ', '', rownames(step.env$anova)) # get variable names from 'ordistep'   and 'envfit' results
  
    Significant <- unique(vars.ordistep)
    df_Significant <- df_Traits[names(df_Traits) %in% vars.ordistep]
    
    #First check of O2 is collinear, otherwise VIF goes up for many values 
    if ("O2" %in% names(which(vif.cca(dbRDA_results) > 10))) {
      print(paste("Removing O2 for Collinearity"))
      df_Significant <- df_Significant[, !names(df_Significant) %in% c("O2")]
     }
    
    # Perform the same RDA but with fewer variables (selected/relevant variables)
    dbRDA_results <- capscale(avgdist_abund_table ~ ., data=df_Significant, distance='bray')
    summary(dbRDA_results, display=NULL) # summary of the results
    anova(dbRDA_results) # permutation test of statistical significance
    
    print("VIF")
    print(sort(vif.cca(dbRDA_results))) # calculate VIF, scores >10 are redundant 
    vars.envfit <- names(which(vif.cca(dbRDA_results) <= 10))

    if (length(names(which(vif.cca(dbRDA_results) > 10))) > 0) {
      print(paste("Removing", names(which(vif.cca(dbRDA_results) > 10)), "for Collinearity"))
    df_Significant <- df_Significant[, !names(df_Significant) %in% names(which(vif.cca(dbRDA_results) > 10))]
    }

    # Perform the same RDA but with fewer variables (selected/relevant variables)
    dbRDA_results <- capscale(avgdist_abund_table ~ ., data=df_Significant, distance='bray')
    summary(dbRDA_results, display=NULL) # summary of the results
    anova(dbRDA_results) # permutation test of statistical significance
    envfit_res <- envfit(dbRDA_results, df_Significant, permutations = 999, na.rm = TRUE)
  
    RDA_list[[RDA_list_length+1]] <- dbRDA_results
    names(RDA_list)[[RDA_list_length+1]] <- paste("dbRDA_results", Datasetname, sep="_")
    
    RDA_list[[RDA_list_length+2]] <- step.env
    names(RDA_list)[[RDA_list_length+2]] <- paste("ordistep_results", Datasetname, sep="_")
    
    RDA_list[[RDA_list_length+3]] <- abund_table.adonis.sig
    names(RDA_list)[[RDA_list_length+3]] <- paste("PERMANOVA_results", Datasetname, sep="_")
    
    RDA_list[[RDA_list_length+4]] <- envfit_res
    names(RDA_list)[[RDA_list_length+4]] <- paste("EnvFit_results", Datasetname, sep="_")
    
    # Set up dataframes for plotting the results
    scores_df <- cbind(df_Traits, scores(dbRDA_results, display="sites"))
    #spp_dpRDA <- cbind(data.frame(tax_table(ps_Filter)), vegan::scores(dbRDA_results, display="species"))
    scores_dbRDA <- vegan::scores(dbRDA_results, display=c("sp","wa","lc","bp","cn"))
  
    ##########
    #Plot RDA#
    ##########
    tryCatch({
    COL <- col.Palette$col.Palette.Season[names(col.Palette$col.Palette.Season) %in% sample_data(ps)$Season]
  
    p <-ggplot(scores_df, aes(x=CAP1, y=CAP2, color=fSeason, shape=fLOC)) + #, shape=RepsSpecs
        ggtitle (paste(Datasetname, "RDA"))+
        scale_color_manual(values = COL) + 
        scale_fill_manual(values = COL) +
        geom_point(size=4.5)

      if (Datasetname %in% c("GC")) {
      p <- p + scale_shape_manual(values=c(16, 17, 8, 25)) 
      } else if (Datasetname %in% c("OE", "WF")) {
      p <- p + scale_shape_manual(values=c(15, 16, 17, 8, 25)) 
      }
      
      #p <- p + ggrepel::geom_text_repel(aes(label = fLOC))
    
      multiplier <- vegan:::ordiArrowMul(scores_dbRDA$biplot)
      df_arrows <- scores_dbRDA$biplot*multiplier
      colnames(df_arrows)<-c("x","y")
      df_arrows=as.data.frame(df_arrows)
    
      p <- p + geom_segment(data=df_arrows,inherit.aes=F, mapping=(aes(x = 0, y = 0, xend = x, yend = y)),
                            arrow = arrow(length = unit(0.2, "cm")),color="grey60",alpha=0.8) +
      geom_text_repel(data=as.data.frame(df_arrows*1.2), inherit.aes=F, mapping=aes(x, y, label =
                      rownames(df_arrows)),color="grey20",alpha=0.99, size = 6)
      
      p <- p + atheme +
      theme(
         panel.background = element_rect(fill='transparent'), #transparent panel bg
         plot.background = element_rect(fill='transparent', color=NA), #transparent plot bg
         ) +
      theme(axis.title.x.bottom =  element_text(color="grey13"), 
        strip.text = element_text(color = "black", face= "bold")) +
      theme(
        legend.title=element_text(size=10), 
        legend.text=element_text(size=10)) +
    theme(
        panel.grid.major = element_line(colour = "white"), 
        panel.grid.minor = element_line(colour = "white"))
  
    # add in Percent to Axis lables 
    p <- p + xlab(paste("RDA1", sub("CAP1", "", ord_labels(dbRDA_results)[1]))) + 
         ylab(paste("RDA2", sub("CAP2", "", ord_labels(dbRDA_results)[2]))) +
      theme(
        axis.title.x.bottom = element_text(size=12,face = "bold"), 
        axis.title.y.left = element_text(size=12,face = "bold"))
    
    # p <- p + annotate("text", x = Inf, y = Inf, hjust = 1, vjust = 1, size = 5,
    #        label = paste(
    #          paste(
    #             "p:",   pslist$PERMANOVA_list[[Dataset]]$Season$Adonis$`Pr(>F)`[1]), 
    #         
    #           paste(
    #            "DF:", pslist$PERMANOVA_list[[Dataset]]$Season$Adonis$Df[3]),
    #          
    #          paste(
    #            "F:",round(pslist$PERMANOVA_list[[Dataset]]$Season$Adonis$F[1])), 
    #          sep = "\n"), 
    #        color = c("darkred", "darkred", "darkred"),  fontface = "bold")

    RDA_plot <- cowplot::plot_grid(p, labels = c(""), ncol = 1)

    ggsave(RDA_plot, filename = paste(paste(save_name, "dbRDA", Datasetname, sep="_"), ".png", sep=""), 
          path = pathPlots, device='png', dpi=300, width = 7,height = 5)
    
    },error=function(e){cat("ERROR :",conditionMessage(e), "\n")})
  
    plot(RDA_plot)
    RDA_list[[RDA_list_length+5]] <- RDA_plot
    names(RDA_list)[[RDA_list_length+5]] <- paste("dbRDA", Datasetname, sep="_")
  }
##  [1] "Salinity"    "Temperature" "Age"         "FCF"         "HSI"        
##  [6] "SSI"         "GSI"         "FillLevel"   "NH4"         "NO2"        
## [11] "NO3"         "O2"          "PO4"         "pH"          "SPM"        
## [16] "TOC"         "fSeason"     "fLOC"       
## [1] "Samples removed for low frequency below 5000 seqs in;ps_OE"
## named integer(0)
## [1] "NAs in TraitData"
## 
## FALSE 
##  2208 
##  [1] Temperature Salinity    FillLevel   Age         HSI         SSI        
##  [7] GSI         FCF         NH4         NO2         NO3         O2         
## [13] PO4         pH          TOC         SPM        
## <0 rows> (or 0-length row.names)
## [1] "Smallest Group for avgdist"
## [1] 9071
## [1] "PERMANOVA OE"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = avgdist_abund_table ~ ., data = df_Include)
##              Df SumOfSqs      R2       F Pr(>F)    
## Temperature   1    3.333 0.10234 25.3357  0.001 ***
## Salinity      1    1.123 0.03448  8.5365  0.001 ***
## FillLevel     1    0.550 0.01689  4.1807  0.002 ** 
## Age           1    0.290 0.00892  2.2070  0.032 *  
## HSI           1    1.185 0.03639  9.0099  0.001 ***
## SSI           1    0.361 0.01108  2.7428  0.013 *  
## GSI           1    3.260 0.10009 24.7785  0.001 ***
## FCF           1    0.597 0.01832  4.5345  0.001 ***
## NH4           1    0.477 0.01463  3.6229  0.007 ** 
## NO2           1    0.911 0.02797  6.9252  0.001 ***
## NO3           1    0.910 0.02794  6.9165  0.001 ***
## O2            1    0.695 0.02133  5.2796  0.001 ***
## PO4           1    1.146 0.03519  8.7124  0.001 ***
## pH            1    0.558 0.01713  4.2398  0.002 ** 
## TOC           1    0.667 0.02047  5.0665  0.001 ***
## SPM           1    0.588 0.01807  4.4725  0.001 ***
## Residual    121   15.920 0.48877                   
## Total       137   32.571 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "VIF"
##         SSI   FillLevel         Age         SPM         HSI         PO4 
##    1.235635    1.238250    1.389237    1.907918    2.037286    2.619247 
##         FCF         TOC         GSI          pH    Salinity         NH4 
##    2.658719    2.989917    3.585188    5.062142    5.551753    7.980814 
##         NO2         NO3 Temperature          O2 
##    8.216231    8.224014   22.591065   27.194257 
## 
## Start: avgdist_abund_table ~ 1 
## 
##               Df    AIC       F Pr(>F)   
## + GSI          1 460.78 28.1692  0.005 **
## + Temperature  1 472.18 15.1623  0.005 **
## + FCF          1 474.40 12.7475  0.005 **
## + SPM          1 475.32 11.7572  0.005 **
## + O2           1 477.93  8.9884  0.005 **
## + NO3          1 479.03  7.8392  0.005 **
## + TOC          1 479.25  7.6061  0.005 **
## + PO4          1 479.68  7.1594  0.005 **
## + HSI          1 479.71  7.1281  0.005 **
## + Salinity     1 481.71  5.0705  0.005 **
## + pH           1 482.40  4.3649  0.005 **
## + NO2          1 483.21  3.5450  0.005 **
## + Age          1 484.02  2.7252  0.010 **
## + NH4          1 484.46  2.2894  0.055 . 
## + FillLevel    1 484.76  1.9892  0.055 . 
## + SSI          1 485.50  1.2535  0.200   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI 
## 
##       Df    AIC      F Pr(>F)   
## - GSI  1 484.76 28.169  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               Df    AIC       F Pr(>F)   
## + Temperature  1 452.19 10.7722  0.005 **
## + O2           1 453.80  9.0860  0.005 **
## + NO3          1 454.29  8.5733  0.005 **
## + FCF          1 456.09  6.7056  0.005 **
## + TOC          1 456.51  6.2787  0.005 **
## + Salinity     1 457.14  5.6404  0.005 **
## + PO4          1 457.21  5.5660  0.005 **
## + SPM          1 457.73  5.0328  0.005 **
## + NO2          1 458.00  4.7602  0.005 **
## + pH           1 458.67  4.0854  0.005 **
## + NH4          1 459.92  2.8358  0.005 **
## + FillLevel    1 459.76  2.9901  0.015 * 
## + HSI          1 460.81  1.9472  0.035 * 
## + Age          1 461.18  1.5752  0.080 . 
## + SSI          1 461.85  0.9143  0.455   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature 
## 
##               Df    AIC      F Pr(>F)   
## - Temperature  1 460.78 10.772  0.005 **
## - GSI          1 472.18 23.315  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             Df    AIC      F Pr(>F)   
## + O2         1 445.64 8.5701  0.005 **
## + PO4        1 448.07 6.0803  0.005 **
## + Salinity   1 448.14 6.0099  0.005 **
## + SPM        1 448.32 5.8262  0.005 **
## + NO2        1 449.13 5.0010  0.005 **
## + pH         1 449.59 4.5392  0.005 **
## + TOC        1 450.14 3.9902  0.005 **
## + NO3        1 450.19 3.9417  0.005 **
## + FCF        1 450.70 3.4284  0.005 **
## + NH4        1 451.03 3.1051  0.005 **
## + FillLevel  1 451.07 3.0616  0.005 **
## + HSI        1 451.86 2.2801  0.015 * 
## + Age        1 452.71 1.4465  0.165   
## + SSI        1 453.72 0.4526  0.915   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 
## 
##               Df    AIC       F Pr(>F)   
## - O2           1 452.19  8.5701  0.005 **
## - Temperature  1 453.80 10.2385  0.005 **
## - GSI          1 463.71 20.9806  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             Df    AIC      F Pr(>F)   
## + pH         1 439.44 8.1328  0.005 **
## + Salinity   1 441.28 6.2722  0.005 **
## + SPM        1 441.30 6.2451  0.005 **
## + TOC        1 442.83 4.7153  0.005 **
## + NO3        1 443.28 4.2677  0.005 **
## + FCF        1 443.70 3.8510  0.005 **
## + PO4        1 443.91 3.6382  0.005 **
## + NO2        1 444.71 2.8483  0.010 **
## + FillLevel  1 445.02 2.5461  0.010 **
## + NH4        1 445.28 2.2890  0.010 **
## + Age        1 446.19 1.4033  0.140   
## + HSI        1 446.50 1.1024  0.290   
## + SSI        1 447.03 0.5882  0.780   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH 
## 
##               Df    AIC       F Pr(>F)   
## - pH           1 445.64  8.1328  0.005 **
## - O2           1 449.59 12.2392  0.005 **
## - Temperature  1 450.23 12.9127  0.005 **
## - GSI          1 453.85 16.7908  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             Df    AIC      F Pr(>F)   
## + SPM        1 434.74 6.5718  0.005 **
## + Salinity   1 436.49 4.8270  0.005 **
## + PO4        1 437.22 4.1023  0.005 **
## + FCF        1 438.11 3.2260  0.005 **
## + NO3        1 438.45 2.8928  0.005 **
## + NH4        1 438.75 2.6028  0.005 **
## + NO2        1 438.76 2.5914  0.005 **
## + TOC        1 438.76 2.5889  0.010 **
## + FillLevel  1 439.57 1.8039  0.035 * 
## + Age        1 439.92 1.4699  0.085 . 
## + HSI        1 440.31 1.0943  0.290   
## + SSI        1 440.77 0.6427  0.800   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM 
## 
##               Df    AIC       F Pr(>F)   
## - SPM          1 439.44  6.5718  0.005 **
## - pH           1 441.30  8.4503  0.005 **
## - O2           1 445.43 12.7130  0.005 **
## - Temperature  1 446.34 13.6720  0.005 **
## - GSI          1 447.19 14.5702  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             Df    AIC      F Pr(>F)   
## + PO4        1 432.19 4.3880  0.005 **
## + Salinity   1 432.31 4.2773  0.005 **
## + NO3        1 433.65 2.9648  0.005 **
## + FCF        1 434.10 2.5304  0.005 **
## + TOC        1 433.78 2.8350  0.010 **
## + NO2        1 434.94 1.7204  0.035 * 
## + NH4        1 434.84 1.8165  0.045 * 
## + FillLevel  1 434.94 1.7228  0.065 . 
## + HSI        1 435.56 1.1214  0.235   
## + Age        1 435.52 1.1629  0.260   
## + SSI        1 436.04 0.6694  0.750   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 
## 
##               Df    AIC       F Pr(>F)   
## - PO4          1 434.74  4.3880  0.010 **
## - SPM          1 437.22  6.8446  0.005 **
## - pH           1 439.27  8.9097  0.005 **
## - O2           1 441.19 10.8699  0.005 **
## - GSI          1 442.52 12.2374  0.005 **
## - Temperature  1 443.00 12.7440  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             Df    AIC      F Pr(>F)   
## + NO3        1 428.65 5.3231  0.005 **
## + Salinity   1 430.59 3.4383  0.005 **
## + TOC        1 430.64 3.3855  0.005 **
## + FCF        1 431.53 2.5314  0.005 **
## + NO2        1 432.35 1.7433  0.020 * 
## + FillLevel  1 432.38 1.7159  0.035 * 
## + NH4        1 432.58 1.5329  0.065 . 
## + Age        1 432.98 1.1488  0.215   
## + HSI        1 433.08 1.0523  0.325   
## + SSI        1 433.56 0.5969  0.835   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 +      NO3 
## 
##               Df    AIC       F Pr(>F)   
## - NO3          1 432.19  5.3231  0.005 **
## - pH           1 432.95  6.0722  0.005 **
## - PO4          1 433.65  6.7608  0.005 **
## - SPM          1 434.12  7.2311  0.005 **
## - GSI          1 434.56  7.6664  0.005 **
## - O2           1 435.75  8.8595  0.005 **
## - Temperature  1 440.64 13.8615  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             Df    AIC      F Pr(>F)   
## + Salinity   1 423.44 6.9246  0.005 **
## + FCF        1 428.22 2.2972  0.005 **
## + TOC        1 428.13 2.3805  0.020 * 
## + NH4        1 428.98 1.5759  0.040 * 
## + NO2        1 428.79 1.7537  0.045 * 
## + FillLevel  1 429.03 1.5267  0.090 . 
## + Age        1 429.38 1.1964  0.225   
## + HSI        1 429.49 1.0972  0.360   
## + SSI        1 429.77 0.8314  0.640   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 +      NO3 + Salinity 
## 
##               Df    AIC       F Pr(>F)   
## - GSI          1 426.94  5.2419  0.005 **
## - SPM          1 427.51  5.7997  0.005 **
## - pH           1 427.78  6.0653  0.005 **
## - Salinity     1 428.65  6.9246  0.005 **
## - PO4          1 429.00  7.2697  0.005 **
## - NO3          1 430.59  8.8445  0.005 **
## - O2           1 431.10  9.3589  0.005 **
## - Temperature  1 437.01 15.4045  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             Df    AIC      F Pr(>F)   
## + TOC        1 419.53 5.5965  0.005 **
## + NH4        1 423.47 1.8371  0.015 * 
## + NO2        1 423.42 1.8848  0.020 * 
## + FCF        1 423.51 1.8049  0.030 * 
## + FillLevel  1 423.77 1.5582  0.120   
## + Age        1 424.13 1.2183  0.210   
## + HSI        1 424.21 1.1473  0.310   
## + SSI        1 424.47 0.8999  0.575   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 +      NO3 + Salinity + TOC 
## 
##               Df    AIC       F Pr(>F)   
## - GSI          1 420.82  3.0869  0.020 * 
## - TOC          1 423.44  5.5965  0.005 **
## - SPM          1 423.74  5.8872  0.005 **
## - PO4          1 424.11  6.2432  0.005 **
## - pH           1 425.51  7.6131  0.005 **
## - NO3          1 425.81  7.9113  0.005 **
## - O2           1 427.14  9.2248  0.005 **
## - Salinity     1 428.13 10.2172  0.005 **
## - Temperature  1 434.06 16.2840  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             Df    AIC      F Pr(>F)  
## + NH4        1 419.22 2.1473  0.015 *
## + NO2        1 419.43 1.9469  0.015 *
## + FillLevel  1 419.83 1.5754  0.060 .
## + Age        1 420.18 1.2535  0.145  
## + HSI        1 420.24 1.1922  0.215  
## + FCF        1 420.41 1.0369  0.325  
## + SSI        1 420.73 0.7391  0.670  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 +      NO3 + Salinity + TOC + NH4 
## 
##               Df    AIC       F Pr(>F)   
## - NH4          1 419.53  2.1473  0.055 . 
## - GSI          1 420.17  2.7486  0.015 * 
## - PO4          1 422.01  4.4859  0.005 **
## - SPM          1 422.29  4.7550  0.005 **
## - TOC          1 423.47  5.8868  0.005 **
## - NO3          1 425.23  7.5899  0.005 **
## - O2           1 425.42  7.7784  0.005 **
## - pH           1 425.57  7.9256  0.005 **
## - Salinity     1 428.08 10.4009  0.005 **
## - Temperature  1 430.99 13.3241  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             Df    AIC      F Pr(>F)   
## + NO2        1 417.65 3.3064  0.005 **
## + FillLevel  1 419.49 1.5879  0.030 * 
## + Age        1 419.85 1.2598  0.170   
## + HSI        1 419.91 1.1990  0.205   
## + FCF        1 420.20 0.9310  0.560   
## + SSI        1 420.41 0.7395  0.685   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 +      NO3 + Salinity + TOC + NH4 + NO2 
## 
##               Df    AIC       F Pr(>F)   
## - GSI          1 418.68  2.7972  0.020 * 
## - NO2          1 419.22  3.3064  0.005 **
## - NH4          1 419.43  3.5074  0.005 **
## - PO4          1 420.70  4.7040  0.005 **
## - SPM          1 420.83  4.8218  0.005 **
## - O2           1 421.78  5.7240  0.005 **
## - TOC          1 422.03  5.9651  0.005 **
## - pH           1 422.73  6.6352  0.005 **
## - NO3          1 423.90  7.7626  0.005 **
## - Salinity     1 426.74 10.5452  0.005 **
## - Temperature  1 428.46 12.2570  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             Df    AIC      F Pr(>F)  
## + FillLevel  1 417.62 1.8441  0.020 *
## + HSI        1 418.07 1.4305  0.100 .
## + Age        1 418.29 1.2315  0.205  
## + FCF        1 418.60 0.9538  0.440  
## + SSI        1 418.71 0.8502  0.590  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 +      NO3 + Salinity + TOC + NH4 + NO2 + FillLevel 
## 
##               Df    AIC       F Pr(>F)   
## - FillLevel    1 417.65  1.8441  0.060 . 
## - GSI          1 418.70  2.8183  0.015 * 
## - NO2          1 419.49  3.5527  0.005 **
## - NH4          1 419.72  3.7627  0.005 **
## - SPM          1 420.46  4.4568  0.005 **
## - PO4          1 420.68  4.6652  0.005 **
## - O2           1 421.09  5.0512  0.005 **
## - TOC          1 422.03  5.9418  0.005 **
## - pH           1 422.11  6.0116  0.005 **
## - NO3          1 423.79  7.6171  0.005 **
## - Salinity     1 426.73 10.4791  0.005 **
## - Temperature  1 427.44 11.1726  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       Df    AIC      F Pr(>F)  
## + HSI  1 417.91 1.5464  0.050 *
## + Age  1 418.19 1.2916  0.140  
## + FCF  1 418.61 0.9150  0.550  
## + SSI  1 418.75 0.7909  0.645  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 +      NO3 + Salinity + TOC + NH4 + NO2 + FillLevel + HSI 
## 
##               Df    AIC       F Pr(>F)   
## - HSI          1 417.62  1.5464  0.125   
## - FillLevel    1 418.07  1.9570  0.055 . 
## - GSI          1 418.72  2.5434  0.015 * 
## - NH4          1 420.22  3.9314  0.010 **
## - NO2          1 420.16  3.8743  0.005 **
## - O2           1 420.77  4.4393  0.005 **
## - SPM          1 420.77  4.4451  0.005 **
## - PO4          1 421.02  4.6751  0.005 **
## - TOC          1 422.37  5.9367  0.005 **
## - pH           1 422.45  6.0105  0.005 **
## - NO3          1 424.24  7.7127  0.005 **
## - Temperature  1 426.68 10.0653  0.005 **
## - Salinity     1 427.16 10.5260  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 +      NO3 + Salinity + TOC + NH4 + NO2 + FillLevel 
## 
##       Df    AIC      F Pr(>F)  
## + HSI  1 417.91 1.5464  0.045 *
## + Age  1 418.19 1.2916  0.215  
## + FCF  1 418.61 0.9150  0.525  
## + SSI  1 418.75 0.7909  0.605  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 +      NO3 + Salinity + TOC + NH4 + NO2 + FillLevel + HSI 
## 
##               Df    AIC       F Pr(>F)   
## - HSI          1 417.62  1.5464  0.105   
## - FillLevel    1 418.07  1.9570  0.065 . 
## - GSI          1 418.72  2.5434  0.025 * 
## - NO2          1 420.16  3.8743  0.010 **
## - NH4          1 420.22  3.9314  0.005 **
## - O2           1 420.77  4.4393  0.005 **
## - SPM          1 420.77  4.4451  0.005 **
## - PO4          1 421.02  4.6751  0.005 **
## - TOC          1 422.37  5.9367  0.005 **
## - pH           1 422.45  6.0105  0.005 **
## - NO3          1 424.24  7.7127  0.005 **
## - Temperature  1 426.68 10.0653  0.005 **
## - Salinity     1 427.16 10.5260  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 +      NO3 + Salinity + TOC + NH4 + NO2 + FillLevel 
## 
##       Df    AIC      F Pr(>F)  
## + HSI  1 417.91 1.5464  0.050 *
## + Age  1 418.19 1.2916  0.125  
## + FCF  1 418.61 0.9150  0.505  
## + SSI  1 418.75 0.7909  0.690  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 +      NO3 + Salinity + TOC + NH4 + NO2 + FillLevel + HSI 
## 
##               Df    AIC       F Pr(>F)   
## - HSI          1 417.62  1.5464  0.145   
## - FillLevel    1 418.07  1.9570  0.075 . 
## - GSI          1 418.72  2.5434  0.035 * 
## - NO2          1 420.16  3.8743  0.005 **
## - NH4          1 420.22  3.9314  0.005 **
## - O2           1 420.77  4.4393  0.005 **
## - SPM          1 420.77  4.4451  0.005 **
## - PO4          1 421.02  4.6751  0.005 **
## - TOC          1 422.37  5.9367  0.005 **
## - pH           1 422.45  6.0105  0.005 **
## - NO3          1 424.24  7.7127  0.005 **
## - Temperature  1 426.68 10.0653  0.005 **
## - Salinity     1 427.16 10.5260  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ GSI + Temperature + O2 + pH + SPM + PO4 +      NO3 + Salinity + TOC + NH4 + NO2 + FillLevel 
## 
##       Df    AIC      F Pr(>F)  
## + HSI  1 417.91 1.5464  0.060 .
## + Age  1 418.19 1.2916  0.125  
## + FCF  1 418.61 0.9150  0.525  
## + SSI  1 418.75 0.7909  0.680  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               Df    AIC       F Pr(>F)   
## - FillLevel    1 417.65  1.8441  0.075 . 
## - GSI          1 418.70  2.8183  0.010 **
## - NO2          1 419.49  3.5527  0.005 **
## - NH4          1 419.72  3.7627  0.005 **
## - SPM          1 420.46  4.4568  0.005 **
## - PO4          1 420.68  4.6652  0.005 **
## - O2           1 421.09  5.0512  0.005 **
## - TOC          1 422.03  5.9418  0.005 **
## - pH           1 422.11  6.0116  0.005 **
## - NO3          1 423.79  7.6171  0.005 **
## - Salinity     1 426.73 10.4791  0.005 **
## - Temperature  1 427.44 11.1726  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       Df    AIC      F Pr(>F)  
## + HSI  1 417.91 1.5464  0.060 .
## + Age  1 418.19 1.2916  0.145  
## + FCF  1 418.61 0.9150  0.580  
## + SSI  1 418.75 0.7909  0.640  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [1] "Ordisep OE"
##               Df    AIC       F Pr(>F)   
## + GSI          1 460.78 28.1692  0.005 **
## + Temperature  1 452.19 10.7722  0.005 **
## + O2           1 445.64  8.5701  0.005 **
## + pH           1 439.44  8.1328  0.005 **
## + SPM          1 434.74  6.5718  0.005 **
## + PO4          1 432.19  4.3880  0.005 **
## + NO3          1 428.65  5.3231  0.005 **
## + Salinity     1 423.44  6.9246  0.005 **
## + TOC          1 419.53  5.5965  0.005 **
## + NH4          1 419.22  2.1473  0.015 * 
## + NO2          1 417.65  3.3064  0.005 **
## + FillLevel    1 417.62  1.8441  0.020 * 
## + HSI          1 417.91  1.5464  0.050 * 
## - HSI          1 417.62  1.5464  0.125   
## + HSI1         1 417.91  1.5464  0.045 * 
## - HSI1         1 417.62  1.5464  0.105   
## + HSI2         1 417.91  1.5464  0.050 * 
## - HSI2         1 417.62  1.5464  0.145   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Removing O2 for Collinearity"
## [1] "VIF"
##   FillLevel         SPM         HSI         PO4         TOC         GSI 
##    1.135031    1.471549    1.778025    2.008939    2.489322    2.839159 
##          pH    Salinity Temperature         NO2         NO3         NH4 
##    3.096458    5.051252    6.326722    6.891286    6.903946    7.898809
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

##  [1] "Salinity"    "Temperature" "Age"         "FCF"         "HSI"        
##  [6] "SSI"         "GSI"         "FillLevel"   "NH4"         "NO2"        
## [11] "NO3"         "O2"          "PO4"         "pH"          "SPM"        
## [16] "TOC"         "fSeason"     "fLOC"       
## [1] "Samples removed for low frequency below 5000 seqs in;ps_GC"
## GCAU21TWFL2 GCAU21MLEB1 
##          45          46 
## [1] "NAs in TraitData"
## 
## FALSE 
##  1376 
##  [1] Temperature Salinity    FillLevel   Age         HSI         SSI        
##  [7] GSI         FCF         NH4         NO2         NO3         O2         
## [13] PO4         pH          TOC         SPM        
## <0 rows> (or 0-length row.names)
## [1] "Smallest Group for avgdist"
## [1] 5998
## [1] "PERMANOVA GC"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = avgdist_abund_table ~ ., data = df_Include)
##             Df SumOfSqs      R2      F Pr(>F)    
## Temperature  1   0.7239 0.04278 4.6766  0.001 ***
## Salinity     1   0.7376 0.04359 4.7650  0.001 ***
## FillLevel    1   0.1426 0.00843 0.9214  0.491    
## Age          1   0.3107 0.01836 2.0070  0.026 *  
## HSI          1   0.1643 0.00971 1.0615  0.339    
## SSI          1   0.0852 0.00504 0.5504  0.908    
## GSI          1   0.1150 0.00680 0.7429  0.719    
## FCF          1   0.1592 0.00941 1.0287  0.377    
## NH4          1   0.5826 0.03444 3.7639  0.001 ***
## NO2          1   1.0029 0.05927 6.4789  0.001 ***
## NO3          1   0.5381 0.03181 3.4766  0.001 ***
## O2           1   0.2650 0.01566 1.7118  0.052 .  
## PO4          1   0.5173 0.03057 3.3418  0.001 ***
## pH           1   0.1817 0.01074 1.1737  0.264    
## TOC          1   0.4240 0.02506 2.7392  0.005 ** 
## SPM          1   0.2886 0.01706 1.8646  0.045 *  
## Residual    69  10.6804 0.63127                  
## Total       85  16.9190 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "VIF"
##         SSI         HSI   FillLevel         FCF         GSI         Age 
##    1.333285    1.637703    1.972504    2.141985    2.198489    2.573010 
##         TOC         SPM         PO4    Salinity         NO3          pH 
##    2.964841    3.298999    5.820725    6.755114    6.913369    7.086578 
##         NO2         NH4 Temperature          O2 
##   13.252686   18.680343   43.279216   48.706396 
## 
## Start: avgdist_abund_table ~ 1 
## 
##               Df    AIC      F Pr(>F)   
## + NO2          1 240.77 6.1761  0.005 **
## + PO4          1 241.76 5.1427  0.005 **
## + O2           1 242.06 4.8424  0.005 **
## + NH4          1 242.83 4.0428  0.005 **
## + Temperature  1 243.14 3.7277  0.005 **
## + Salinity     1 243.25 3.6144  0.005 **
## + SPM          1 243.32 3.5498  0.005 **
## + TOC          1 244.10 2.7520  0.010 **
## + Age          1 244.36 2.4953  0.010 **
## + pH           1 244.09 2.7618  0.015 * 
## + NO3          1 244.48 2.3697  0.015 * 
## + GSI          1 244.95 1.9058  0.055 . 
## + HSI          1 245.71 1.1420  0.290   
## + FillLevel    1 246.01 0.8504  0.575   
## + SSI          1 245.96 0.8968  0.580   
## + FCF          1 246.12 0.7377  0.720   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ NO2 
## 
##       Df    AIC      F Pr(>F)   
## - NO2  1 244.88 6.1761  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               Df    AIC      F Pr(>F)   
## + PO4          1 236.46 6.3245  0.005 **
## + O2           1 239.22 3.5021  0.005 **
## + Temperature  1 239.22 3.5001  0.005 **
## + TOC          1 239.50 3.2160  0.005 **
## + NO3          1 240.00 2.7179  0.005 **
## + SPM          1 240.06 2.6559  0.005 **
## + Salinity     1 240.26 2.4608  0.015 * 
## + pH           1 240.25 2.4675  0.020 * 
## + NH4          1 240.35 2.3687  0.020 * 
## + Age          1 240.50 2.2198  0.020 * 
## + GSI          1 241.10 1.6293  0.095 . 
## + HSI          1 241.53 1.2127  0.220   
## + FillLevel    1 241.70 1.0425  0.415   
## + FCF          1 242.01 0.7424  0.700   
## + SSI          1 242.04 0.7132  0.705   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ NO2 + PO4 
## 
##       Df    AIC      F Pr(>F)   
## - PO4  1 240.77 6.3245  0.005 **
## - NO2  1 241.76 7.3601  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               Df    AIC      F Pr(>F)   
## + TOC          1 235.16 3.2023  0.005 **
## + Temperature  1 235.23 3.1358  0.005 **
## + pH           1 235.55 2.8241  0.010 **
## + O2           1 235.62 2.7539  0.010 **
## + NO3          1 235.64 2.7311  0.010 **
## + Salinity     1 235.86 2.5150  0.015 * 
## + Age          1 235.99 2.3840  0.015 * 
## + SPM          1 236.24 2.1419  0.015 * 
## + GSI          1 236.78 1.6143  0.035 * 
## + NH4          1 236.84 1.5590  0.085 . 
## + FillLevel    1 237.10 1.3021  0.120   
## + HSI          1 237.22 1.1887  0.260   
## + SSI          1 237.72 0.7051  0.715   
## + FCF          1 237.71 0.7135  0.810   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ NO2 + PO4 + TOC 
## 
##       Df    AIC      F Pr(>F)   
## - TOC  1 236.46 3.2023  0.005 **
## - PO4  1 239.50 6.2742  0.005 **
## - NO2  1 241.71 8.5623  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               Df    AIC      F Pr(>F)   
## + Salinity     1 233.51 3.5185  0.005 **
## + pH           1 233.68 3.3516  0.005 **
## + O2           1 234.22 2.8225  0.005 **
## + Age          1 234.55 2.5006  0.005 **
## + NO3          1 234.70 2.3561  0.005 **
## + Temperature  1 233.85 3.1855  0.010 **
## + SPM          1 235.12 1.9481  0.035 * 
## + GSI          1 235.37 1.7065  0.075 . 
## + HSI          1 235.87 1.2312  0.215   
## + NH4          1 235.84 1.2544  0.235   
## + FillLevel    1 235.94 1.1571  0.295   
## + FCF          1 236.37 0.7551  0.730   
## + SSI          1 236.52 0.6071  0.850   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ NO2 + PO4 + TOC + Salinity 
## 
##            Df    AIC      F Pr(>F)   
## - TOC       1 235.86 4.2059  0.010 **
## - Salinity  1 235.16 3.5185  0.005 **
## - PO4       1 238.44 6.8023  0.005 **
## - NO2       1 238.46 6.8193  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               Df    AIC      F Pr(>F)   
## + Temperature  1 231.72 3.6029  0.005 **
## + O2           1 232.39 2.9503  0.005 **
## + pH           1 232.59 2.7648  0.005 **
## + NO3          1 233.07 2.3034  0.010 **
## + SPM          1 233.43 1.9550  0.010 **
## + Age          1 233.03 2.3350  0.015 * 
## + GSI          1 233.57 1.8222  0.060 . 
## + NH4          1 234.13 1.2891  0.135   
## + HSI          1 234.12 1.2965  0.175   
## + FillLevel    1 234.15 1.2734  0.210   
## + FCF          1 234.77 0.6894  0.860   
## + SSI          1 234.86 0.6041  0.865   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ NO2 + PO4 + TOC + Salinity + Temperature 
## 
##               Df    AIC      F Pr(>F)   
## - Temperature  1 233.51 3.6029  0.005 **
## - Salinity     1 233.85 3.9336  0.005 **
## - TOC          1 234.41 4.4878  0.005 **
## - NO2          1 236.18 6.2430  0.005 **
## - PO4          1 236.52 6.5854  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             Df    AIC      F Pr(>F)   
## + NO3        1 230.67 2.8494  0.010 **
## + Age        1 231.25 2.2975  0.015 * 
## + NH4        1 231.79 1.7909  0.025 * 
## + O2         1 231.94 1.6460  0.055 . 
## + SPM        1 232.26 1.3523  0.160   
## + pH         1 232.47 1.1578  0.260   
## + FillLevel  1 232.61 1.0224  0.380   
## + HSI        1 232.62 1.0152  0.445   
## + FCF        1 232.85 0.7977  0.695   
## + GSI        1 233.00 0.6620  0.820   
## + SSI        1 233.06 0.6115  0.900   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ NO2 + PO4 + TOC + Salinity + Temperature +      NO3 
## 
##               Df    AIC      F Pr(>F)   
## - NO3          1 231.72 2.8494  0.005 **
## - TOC          1 232.48 3.5741  0.005 **
## - Salinity     1 232.69 3.7767  0.005 **
## - PO4          1 232.87 3.9555  0.005 **
## - Temperature  1 233.07 4.1417  0.005 **
## - NO2          1 234.53 5.5672  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             Df    AIC      F Pr(>F)  
## + Age        1 230.67 1.8337  0.035 *
## + SPM        1 231.17 1.3753  0.080 .
## + NH4        1 230.97 1.5560  0.115  
## + O2         1 231.19 1.3597  0.170  
## + pH         1 231.41 1.1531  0.275  
## + HSI        1 231.48 1.0886  0.365  
## + FillLevel  1 231.71 0.8739  0.575  
## + FCF        1 231.84 0.7592  0.635  
## + GSI        1 231.93 0.6804  0.855  
## + SSI        1 232.03 0.5823  0.930  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ NO2 + PO4 + TOC + Salinity + Temperature +      NO3 + Age 
## 
##               Df    AIC      F Pr(>F)   
## - Age          1 230.67 1.8337  0.060 . 
## - NO3          1 231.25 2.3756  0.015 * 
## - Salinity     1 231.62 2.7195  0.015 * 
## - TOC          1 231.69 2.7842  0.015 * 
## - PO4          1 232.29 3.3534  0.005 **
## - Temperature  1 232.55 3.5947  0.005 **
## - NO2          1 234.17 5.1462  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##             Df    AIC      F Pr(>F)
## + NH4        1 231.06 1.4530  0.105
## + O2         1 231.11 1.4141  0.125
## + SPM        1 231.12 1.4041  0.185
## + HSI        1 231.44 1.1110  0.230
## + pH         1 231.38 1.1634  0.290
## + FillLevel  1 231.93 0.6668  0.840
## + SSI        1 231.99 0.6167  0.865
## + FCF        1 232.04 0.5713  0.920
## + GSI        1 232.08 0.5327  0.980
## 
## [1] "Ordisep GC"
##               Df    AIC      F Pr(>F)   
## + NO2          1 240.77 6.1761  0.005 **
## + PO4          1 236.46 6.3245  0.005 **
## + TOC          1 235.16 3.2023  0.005 **
## + Salinity     1 233.51 3.5185  0.005 **
## + Temperature  1 231.72 3.6029  0.005 **
## + NO3          1 230.67 2.8494  0.010 **
## + Age          1 230.67 1.8337  0.035 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Removing O2 for Collinearity"
## [1] "VIF"
##         NO2         Age         TOC    Salinity         PO4         NO3 
##    1.454990    1.725781    1.774548    1.851783    2.344520    4.051192 
## Temperature 
##    4.052744
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

##  [1] "Salinity"    "Temperature" "NH4"         "NO2"         "NO3"        
##  [6] "O2"          "PO4"         "pH"          "SPM"         "TOC"        
## [11] "fSeason"     "fLOC"       
## [1] "No Waterfiler Samples removed ;ps_WF"
## [1] "NAs in TraitData"
## 
## FALSE 
##   190 
##  [1] Temperature Salinity    NH4         NO2         NO3         O2         
##  [7] PO4         pH          TOC         SPM        
## <0 rows> (or 0-length row.names)
## [1] "Smallest Group for avgdist"
## [1] 3576
## [1] "PERMANOVA WF"
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## vegan::adonis2(formula = avgdist_abund_table ~ ., data = df_Include)
##             Df SumOfSqs      R2      F Pr(>F)    
## Temperature  1   0.8889 0.20631 8.5455  0.001 ***
## Salinity     1   0.8819 0.20467 8.4778  0.001 ***
## NH4          1   0.2410 0.05593 2.3166  0.028 *  
## NO2          1   0.1149 0.02667 1.1046  0.366    
## NO3          1   0.2782 0.06456 2.6741  0.011 *  
## O2           1   0.2057 0.04774 1.9775  0.042 *  
## PO4          1   0.2274 0.05279 2.1865  0.029 *  
## pH           1   0.3006 0.06977 2.8898  0.013 *  
## TOC          1   0.1443 0.03348 1.3869  0.187    
## SPM          1   0.1937 0.04495 1.8620  0.085 .  
## Residual     8   0.8322 0.19314                  
## Total       18   4.3087 1.00000                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "VIF"
##    Salinity         PO4          pH         TOC         NO3         SPM 
##    7.176865    9.171618    9.915326   24.715622   47.988849   53.281803 
##         NO2         NH4 Temperature          O2 
##   58.507287  108.432044  254.059907  278.295693 
## 
## Start: avgdist_abund_table ~ 1 
## 
##               Df    AIC      F Pr(>F)   
## + pH           1 25.785 5.0481  0.005 **
## + Salinity     1 26.306 4.4512  0.005 **
## + Temperature  1 26.335 4.4188  0.005 **
## + NO3          1 26.665 4.0497  0.005 **
## + O2           1 27.159 3.5098  0.005 **
## + TOC          1 28.332 2.2816  0.025 * 
## + SPM          1 28.572 2.0399  0.045 * 
## + PO4          1 29.051 1.5652  0.095 . 
## + NO2          1 29.343 1.2821  0.175   
## + NH4          1 29.669 0.9711  0.440   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ pH 
## 
##      Df    AIC      F Pr(>F)   
## - pH  1 28.725 5.0481  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               Df    AIC      F Pr(>F)   
## + Temperature  1 21.121 6.7209  0.005 **
## + O2           1 21.776 5.9512  0.005 **
## + NO3          1 22.398 5.2446  0.005 **
## + TOC          1 24.965 2.5599  0.015 * 
## + Salinity     1 25.726 1.8311  0.055 . 
## + PO4          1 25.908 1.6608  0.090 . 
## + SPM          1 26.143 1.4436  0.120   
## + NH4          1 26.759 0.8874  0.510   
## + NO2          1 27.094 0.5923  0.940   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ pH + Temperature 
## 
##               Df    AIC      F Pr(>F)   
## - Temperature  1 25.785 6.7209  0.005 **
## - pH           1 26.335 7.3885  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            Df    AIC      F Pr(>F)   
## + PO4       1 20.503 2.2162  0.005 **
## + NO3       1 20.669 2.0664  0.010 **
## + Salinity  1 21.280 1.5260  0.055 . 
## + TOC       1 21.752 1.1207  0.270   
## + NH4       1 21.933 0.9685  0.560   
## + SPM       1 22.003 0.9096  0.565   
## + NO2       1 22.135 0.7989  0.755   
## + O2        1 22.275 0.6835  0.905   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ pH + Temperature + PO4 
## 
##               Df    AIC      F Pr(>F)   
## - PO4          1 21.121 2.2162  0.030 * 
## - Temperature  1 25.908 7.1490  0.005 **
## - pH           1 26.208 7.5014  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            Df    AIC      F Pr(>F)  
## + Salinity  1 20.131 1.8615  0.015 *
## + TOC       1 20.584 1.4878  0.070 .
## + NO3       1 20.690 1.4019  0.145  
## + SPM       1 21.001 1.1515  0.310  
## + NO2       1 21.152 1.0314  0.400  
## + NH4       1 21.364 0.8649  0.620  
## + O2        1 21.580 0.6970  0.860  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: avgdist_abund_table ~ pH + Temperature + PO4 + Salinity 
## 
##               Df    AIC      F Pr(>F)   
## - Salinity     1 20.503 1.8615  0.065 . 
## - PO4          1 21.280 2.5240  0.010 **
## - pH           1 22.071 3.2262  0.005 **
## - Temperature  1 25.327 6.4457  0.005 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       Df    AIC      F Pr(>F)  
## + NO3  1 20.052 1.5030  0.060 .
## + NO2  1 20.569 1.1143  0.300  
## + SPM  1 20.556 1.1233  0.305  
## + TOC  1 20.637 1.0637  0.380  
## + NH4  1 20.753 0.9780  0.425  
## + O2   1 21.073 0.7443  0.775  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## [1] "Ordisep WF"
##               Df    AIC      F Pr(>F)   
## + pH           1 25.785 5.0481  0.005 **
## + Temperature  1 21.121 6.7209  0.005 **
## + PO4          1 20.503 2.2162  0.005 **
## + Salinity     1 20.131 1.8615  0.015 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Removing O2 for Collinearity"
## [1] "VIF"
## Temperature         PO4          pH    Salinity 
##    1.197473    1.199426    3.014100    3.091885
## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

pslist$RDA_list <- RDA_list
print(sort(vif.cca(RDA_list$dbRDA_results_OE)))
##   FillLevel         SPM         HSI         PO4         TOC         GSI 
##    1.135031    1.471549    1.778025    2.008939    2.489322    2.839159 
##          pH    Salinity Temperature         NO2         NO3         NH4 
##    3.096458    5.051252    6.326722    6.891286    6.903946    7.898809
print(sort(vif.cca(RDA_list$dbRDA_results_GC)))
##         NO2         Age         TOC    Salinity         PO4         NO3 
##    1.454990    1.725781    1.774548    1.851783    2.344520    4.051192 
## Temperature 
##    4.052744
#Export EnvFit Results
for (Dataset in c(names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE", names(pslist)[!grepl("WF", names(pslist))])], "ps_WF")) {
  Datasetname <- sub("ps_", "", Dataset)
  
  EnvFit_results<- as.data.frame(
    cbind(round(pslist$RDA_list[[paste("EnvFit_results", Datasetname, sep="_")]]$vectors$r, digits=3), 
    round(pslist$RDA_list[[paste("EnvFit_results", Datasetname, sep="_")]]$vectors$pvals, digits=3))
    ) %>% 
    setNames(c("R2", "P"))

  write.csv2(EnvFit_results, paste0(file.path(path_Output_16S, paste(paste(save_name, 
            paste("EnvFit_results", sep=""),  Datasetname, sep="_"),".csv", sep=""))))
}

6.3.5 Mantel Test

https://jkzorz.github.io/2019/07/08/mantel-test.html Mantel tests are correlation tests that determine the correlation between two matrices (rather than two variables). When using the test for microbial ecology, the matrices are often distance/dissimilarity matrices with corresponding positions (i.e. samples in the same order in both matrices). In order to calculate the correlation, the matrix values of both matrices are ‘unfolded’ into long column vectors, which are then used to determine correlation. Permutations of one matrix are used to determine significance. Learn more about Mantel tests here.

#########################   
#Fish_OTU vs Environment#
#########################
#based in the avgdist matrix
Mantel_list <- list()
  for (Dataset in names(pslist$BrayCurtis)[grepl("OE_Avg_dist_matrix|GC_Avg_dist_matrix", names(pslist$BrayCurtis ))]) {
  
  Mantel_list_length <- length(Mantel_list)
  Datasetname      <- sub("_Avg_dist_matrix", "", Dataset)
  ps <- pslist[[paste("ps", Datasetname, sep="_")]]

  ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
  
  metadata <- as.data.frame(sample_data(ps_Filter))
  metadata <- metadata[,names(metadata) %in% c("SampleID", "Temperature", "Salinity", "O2", "SPM", "PO4", "NH4", "NO3", "NO2")]


  dist_tibble <- pslist$BrayCurtis[[Dataset]] %>%
    as.matrix() %>%
    as_tibble(rownames="sample") %>%
    pivot_longer(-sample)

  mantel_function = function(species, parameter){
  
  rare_dist_matrix <- dist_tibble %>% 
    left_join(., metadata, by = c("sample" = "SampleID"))%>% 
    left_join(., metadata, by = c("name" = "SampleID")) %>% 
    select(sample, name, value) %>%
    pivot_wider(names_from = name, values_from = value) %>% 
    select(-sample) %>% 
    as.dist()  
  
  env_dist <- dist_tibble %>% 
    left_join(., metadata, by = c("sample" = "SampleID")) %>% 
    left_join(., metadata, by = c("name" = "SampleID")) %>% 
    select(sample, paste0(parameter,".x")) %>%
    unique() %>%
    column_to_rownames("sample") %>%
    dist(method="euclidean")
  
  abund_temp = mantel(rare_dist_matrix, env_dist, method = "spearman", permutations = 9999, na.rm = TRUE)
  cbind(species=species, env=parameter, mantel=abund_temp$statistic, sig=abund_temp$signif)
}

  
  
  Env.out = data.frame(species=NA, env=NA, mantel=NA, sig=NA)
  species=Datasetname
  envs=c("Temperature", "Salinity", "O2", "SPM", "PO4", "NH4", "NO3", "NO2")
  for(Dat in species){
    for(params in envs)
      Env.out <- rbind(Env.out, mantel_function(species = Dat, parameter = params))
  }
  Mantel_list[[Mantel_list_length+1]] <- Env.out
  names(Mantel_list)[[Mantel_list_length+1]] <- paste(Datasetname, "Env", sep="_")
}  

########################   
#Fish_OTU vs Physiology#
########################
#based in the avgdist matrix
for (Dataset in names(pslist$BrayCurtis)[grepl("OE_Avg_dist_matrix|GC_Avg_dist_matrix", names(pslist$BrayCurtis))]) {
  
  Mantel_list_length <- length(Mantel_list)
  Datasetname      <- sub("_Avg_dist_matrix", "", Dataset)
  ps <- pslist[[paste("ps", Datasetname, sep="_")]]

  ps_Filter <- subset_samples(ps, sample_sums(ps) > 5000)
  
  Phys <- c(
   "Length", "Age", "FCF", "GSI", "HSI", "SSI", "FillLevel")

  '%!in%' <- function(x,y)!('%in%'(x,y))
  
  df_Traits <- data.frame(sample_data(ps_Filter))
  df_Traits <- df_Traits[names(df_Traits) %in% Phys]

  print("NAs in TraitData")
  print(table(is.na(df_Traits)))
  
  print((df_Traits[!complete.cases(df_Traits), ]))
  
  df_Traits  <- df_Traits[complete.cases(df_Traits), ]
  
  df_Traits <- decostand(df_Traits, method='standardize')
  
  ps_Filter <- prune_samples(sample_names(ps_Filter)[sample_names(ps_Filter) %in%
               rownames(df_Traits[complete.cases(df_Traits), ])], ps_Filter)

  ##################################################
  #Determine which variables are useful for the plot#
  ##################################################
  abund_table <- as.data.frame(otu_table(ps_Filter, taxa_are_rows = FALSE)) 
  smallest_group <- min(rowSums(otu_table(ps_Filter)))
  avgdist_abund_table  <- avgdist(abund_table, dmethod = "bray", sample = smallest_group)
  
  dist_tibble <- avgdist_abund_table %>%
    as.matrix() %>%
    as_tibble(rownames="sample") %>%
    pivot_longer(-sample)
  
  metadata <- as.data.frame(sample_data(ps_Filter))
  metadata <- metadata[,names(metadata) %in% c("SampleID", Phys)]
  
  mantel_function = function(species, parameter){
  
  rare_dist_matrix <- dist_tibble %>% 
    left_join(., metadata, by = c("sample" = "SampleID"))%>% 
    left_join(., metadata, by = c("name" = "SampleID")) %>% 
    select(sample, name, value) %>%
    pivot_wider(names_from = name, values_from = value) %>% 
    select(-sample) %>% 
    as.dist()  
  
  env_dist <- dist_tibble %>% 
    left_join(., metadata, by = c("sample" = "SampleID")) %>% 
    left_join(., metadata, by = c("name" = "SampleID")) %>% 
    select(sample, paste0(parameter,".x")) %>%
    unique() %>%
    column_to_rownames("sample") %>%
    dist(method="euclidean")
  
  abund_temp = mantel(rare_dist_matrix, env_dist, method = "spearman", permutations = 9999, na.rm = TRUE)
  cbind(species=species, env=parameter, mantel=abund_temp$statistic, sig=abund_temp$signif)
  }

  Phys.out = data.frame(species=NA, env=NA, mantel=NA, sig=NA)
  species=Datasetname
  for(Dat in species){
    for(params in Phys)
      Phys.out <- rbind(Phys.out, mantel_function(species = Dat, parameter = params))
  }
  
  Mantel_list[[Mantel_list_length+1]] <- Phys.out
  names(Mantel_list)[[Mantel_list_length+1]] <- paste(Datasetname, "Phys", sep="_")
}  
## [1] "NAs in TraitData"
## 
## FALSE 
##   966 
## [1] Length    FillLevel Age       HSI       SSI       GSI       FCF      
## <0 rows> (or 0-length row.names)
## [1] "NAs in TraitData"
## 
## FALSE 
##   602 
## [1] Length    FillLevel Age       HSI       SSI       GSI       FCF      
## <0 rows> (or 0-length row.names)
#########################   
#Water_OTU vs Environment#
#########################
#based in the avgdist matrix
for (Dataset in names(pslist)[grepl("ps_WF", names(pslist))]) {
  
  Mantel_list_length <- length(Mantel_list)
  Datasetname      <- sub("ps_", "", Dataset)
  
  ps <- pslist[[paste("ps", Datasetname, sep="_")]]
  
  ps_Filter <- ps
  
  print(paste("Samples removed for low frequency below 5000 seqs in;", Dataset, sep = ""))
  print(which(rowSums(otu_table(ps)) < 5000))
  print("Waterfilters will not be filtered")

  Envs <- c(
   "Temperature", "Salinity", "O2", "SPM", "PO4", "NH4", "NO3", "NO2")


  '%!in%' <- function(x,y)!('%in%'(x,y))
  
  df_Traits <- data.frame(sample_data(ps_Filter))
  df_Traits <- df_Traits[names(df_Traits) %in% Envs]

  print("NAs in TraitData")
  print(table(is.na(df_Traits)))
  
  print((df_Traits[!complete.cases(df_Traits), ]))
  
  df_Traits  <- df_Traits[complete.cases(df_Traits), ]
  
  df_Traits <- decostand(df_Traits, method='standardize')
  
  ps_Filter <- prune_samples(sample_names(ps_Filter)[sample_names(ps_Filter) %in%
               rownames(df_Traits[complete.cases(df_Traits), ])], ps_Filter)

  ##################################################
  #Determine which variables are useful for the plot#
  ##################################################
  abund_table <- as.data.frame(otu_table(ps_Filter, taxa_are_rows = FALSE)) 
  smallest_group <- min(rowSums(otu_table(ps_Filter)))
  print("Smallest group for avgdist")
  print(smallest_group)
  avgdist_abund_table  <- avgdist(abund_table, dmethod = "bray", sample = smallest_group)
  
  dist_tibble <- avgdist_abund_table %>%
    as.matrix() %>%
    as_tibble(rownames="sample") %>%
    pivot_longer(-sample)
  
  metadata <- as.data.frame(sample_data(ps_Filter))
  metadata <- metadata[,names(metadata) %in% c("SampleID", Envs)]

  mantel_function = function(species, parameter){
  
  rare_dist_matrix <- dist_tibble %>% 
    left_join(., metadata, by = c("sample" = "SampleID"))%>% 
    left_join(., metadata, by = c("name" = "SampleID")) %>% 
    select(sample, name, value) %>%
    pivot_wider(names_from = name, values_from = value) %>% 
    select(-sample) %>% 
    as.dist()  
  
  env_dist <- dist_tibble %>% 
    left_join(., metadata, by = c("sample" = "SampleID")) %>% 
    left_join(., metadata, by = c("name" = "SampleID")) %>% 
    select(sample, paste0(parameter,".x")) %>%
    unique() %>%
    column_to_rownames("sample") %>%
    dist(method="euclidean")
  
  abund_temp = mantel(rare_dist_matrix, env_dist, method = "spearman", permutations = 9999, na.rm = TRUE)
  cbind(species=species, env=parameter, mantel=abund_temp$statistic, sig=abund_temp$signif)
}

  Env.out = data.frame(species=NA, env=NA, mantel=NA, sig=NA)
  species=Datasetname
  envs=Envs
  for(Dat in species){
    for(params in envs)
      Env.out <- rbind(Env.out, mantel_function(species = Dat, parameter = params))
  }
  Mantel_list[[Mantel_list_length+1]] <- Env.out
  names(Mantel_list)[[Mantel_list_length+1]] <- paste(Datasetname, "Env", sep="_")
}  
## [1] "Samples removed for low frequency below 5000 seqs in;ps_WF"
## WFSU21SSFL WFWI22MGEB 
##          3          6 
## [1] "Waterfilters will not be filtered"
## [1] "NAs in TraitData"
## 
## FALSE 
##   152 
## [1] Temperature Salinity    NH4         NO2         NO3         O2         
## [7] PO4         SPM        
## <0 rows> (or 0-length row.names)
## [1] "Smallest group for avgdist"
## [1] 3576
#######################
#Fish_ASV vs Water_ASV#
#######################
#-> Problem No Waterfilters for Autumn and SPTW
#Distance Matrix Fish 
#Distance Matrix Water
  #->Elongate Water Matrix to fit the fish Matrix

# for (Dataset in names(pslist$BrayCurtis)[grepl("OE_Avg_dist_matrix|GC_Avg_dist_matrix", names(pslist$BrayCurtis ))]) {
#   
#   Mantel_list_length <- length(Mantel_list)
#   Datasetname      <- sub("_Avg_dist_matrix", "", Dataset)
#   
#   ps_Fish <- pslist[[paste("ps", Datasetname, sep="_")]]
#   ps_Filter_Fish <- subset_samples(ps_Fish, sample_sums(ps_Fish) > 5000)
#   
#   ps_WF <- pslist[[paste("ps", "WF", sep="_")]]
#   abund_table_WF <- data.frame(otu_table(ps_WF))
# 
#   ###########################
#   #Avg_dist for Fish samples#
#   ###########################
#   abund_table_Fish <- as.data.frame(otu_table(ps_Filter_Fish, taxa_are_rows = FALSE)) 
#   print(Datasetname)
#   print("Smalles Group AVG dist")
#   smallest_group_Fish <- min(rowSums(otu_table(ps_Filter_Fish)))
#   print(smallest_group_Fish)
#   avgdist_abund_table_Fish  <- avgdist(abund_table_Fish, dmethod = "bray", sample = smallest_group_Fish)
#   
#   dist_tibble_Fish <- avgdist_abund_table_Fish %>%
#     as.matrix() %>%
#     as_tibble(rownames="sample") %>%
#     pivot_longer(-sample)
# 
#   ################################################
#   #Avg_fist for artificially matched Watersamples#
#   ################################################
#   Names_Fish <- substring(sample_names(ps_Filter_Fish), first = 3)
#   WF_mantel<- list()
#   for (i in seq_along(Names_Fish)) {
#         NAME  <- Names_Fish[i]
#         NAME2 <- gsub("\\d+$", "", NAME)
#         NAME3 <- substr(NAME2 , 1, nchar(NAME2 ) - 2)
#         NAME4 <- substr(NAME3 , 1, nchar(NAME3 ) - 2)
# 
#          if (paste("WF", NAME3, "EB", sep="") %in% rownames(abund_table_WF)) {
#             WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME3, "EB",
#                                                                                  sep=""),]
#           
#         } else if (paste("WF", NAME3, "FL", sep="") %in% rownames(abund_table_WF)){
#             WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME3, "FL",
#                                                                                  sep=""),]
#           
#         } else if (paste("WF", NAME4, "TWEB", sep="") %in% rownames(abund_table_WF)){
#             WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME4, "TWEB",
#                                                                                  sep=""),]
#           
#         } else if (paste("WF", NAME4, "MLFL", sep="") %in% rownames(abund_table_WF)){
#             WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME4, "MLFL",
#                                                                                  sep=""),]
#           
#         } else {print(NA)}
#         
#         names(WF_mantel)[[i]] <- paste("WF", NAME, sep="")
#     }
# 
#   WF_mantel_df <- do.call(rbind, WF_mantel)  
# 
#   print("Smallest Group AVG dist Waterfilter")
#   smallest_group_WF <- min(rowSums(otu_table(ps_WF)))
#   print(smallest_group_WF)
#   avgdist_abund_table_WF  <- avgdist(WF_mantel_df , dmethod = "bray", sample = smallest_group_WF)
#   
#   #####################
#   #Run the Mantel Test#
#   #####################
#   dist_tibble_Fish <- avgdist_abund_table_Fish %>%
#     as.matrix() %>%
#     as_tibble(rownames="sample") %>%
#     pivot_longer(-sample)
#   
#   dist_tibble_WF <- avgdist_abund_table_WF %>%
#     as.matrix() %>%
#     as_tibble(rownames="sample") %>%
#     pivot_longer(-sample)
# 
#   rare_dist_matrix_Fish <- dist_tibble_Fish %>% 
#     left_join(., metadata, by = c("sample" = "SampleID"))%>% 
#     left_join(., metadata, by = c("name" = "SampleID")) %>% 
#     select(sample, name, value) %>%
#     pivot_wider(names_from = name, values_from = value) %>% 
#     select(-sample) %>% 
#     as.dist()  
#   
#   rare_dist_matrix_WF <- dist_tibble_WF %>% 
#     left_join(., metadata, by = c("sample" = "SampleID"))%>% 
#     left_join(., metadata, by = c("name" = "SampleID")) %>% 
#     select(sample, name, value) %>%
#     pivot_wider(names_from = name, values_from = value) %>% 
#     select(-sample) %>% 
#     as.dist()  
# 
#   abund_temp <- 
#     mantel(rare_dist_matrix_Fish, rare_dist_matrix_WF, method = "spearman", permutations = 9999, na.rm = TRUE)
#   Env.out = data.frame(species=Datasetname, env="WF", mantel=abund_temp$statistic, sig=abund_temp$signif)
#   
#     Mantel_list[[Mantel_list_length+1]] <- Env.out
#     names(Mantel_list)[[Mantel_list_length+1]] <- paste(Datasetname, "WF", sep="_")
# }  


##########################################################################
#Without Autumn Samples as there is a large lack in matching watersamples#
##########################################################################
for (Dataset in names(pslist$BrayCurtis)[grepl("OE_Avg_dist_matrix|GC_Avg_dist_matrix", names(pslist$BrayCurtis ))]) {

  Mantel_list_length <- length(Mantel_list)
  Datasetname      <- sub("_Avg_dist_matrix", "", Dataset)

  ps_Fish <- pslist[[paste("ps", Datasetname, sep="_")]]
  ps_Filter_Fish <- subset_samples(ps_Fish, sample_sums(ps_Fish) > 5000)

  ps_Filter_Fish <- prune_samples(!grepl("AU21BB|AU21MG|AU21ML|AU21SS|SU21TW|WI22TW|SP22TW", sample_names(ps_Filter_Fish)), ps_Filter_Fish)


  ps_WF <- pslist[[paste("ps", "WF", sep="_")]]
  abund_table_WF <- data.frame(otu_table(ps_WF))

  ###########################
  #Avg_dist for Fish samples#
  ###########################
  abund_table_Fish <- as.data.frame(otu_table(ps_Filter_Fish, taxa_are_rows = FALSE))
  print(Datasetname)
  print("Smalles Group AVG dist")
  smallest_group_Fish <- min(rowSums(otu_table(ps_Filter_Fish)))
  print(smallest_group_Fish)
  avgdist_abund_table_Fish  <- avgdist(abund_table_Fish, dmethod = "bray", sample = smallest_group_Fish)

  dist_tibble_Fish <- avgdist_abund_table_Fish %>%
    as.matrix() %>%
    as_tibble(rownames="sample") %>%
    pivot_longer(-sample)

  ################################################
  #Avg_fist for artificially matched Watersamples#
  ################################################
  #we lack some watersamples matching with all fish data, as the bacterioplankton appears rather homogenous compared to the fish, we artificially match water samples from the closest time and space with the fish data
  Names_Fish <- substring(sample_names(ps_Filter_Fish), first = 3)
  WF_mantel<- list()
  for (i in seq_along(Names_Fish)) {
        NAME  <- Names_Fish[i]
        NAME2 <- gsub("\\d+$", "", NAME)
        NAME3 <- substr(NAME2 , 1, nchar(NAME2 ) - 2)
        NAME4 <- substr(NAME3 , 1, nchar(NAME3 ) - 2)

         if (paste("WF", NAME3, "EB", sep="") %in% rownames(abund_table_WF)) {
            WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME3, "EB",
                                                                                 sep=""),]

        } else if (paste("WF", NAME3, "FL", sep="") %in% rownames(abund_table_WF)){
            WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME3, "FL",
                                                                                 sep=""),]

        } else if (paste("WF", NAME4, "TWEB", sep="") %in% rownames(abund_table_WF)){
            WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME4, "TWEB",
                                                                                 sep=""),]

        } else if (paste("WF", NAME4, "MLFL", sep="") %in% rownames(abund_table_WF)){
            WF_mantel[[i]] <- abund_table_WF[rownames(abund_table_WF) %in% paste("WF", NAME4, "MLFL",
                                                                                 sep=""),]

        } else {print(NA)}

        names(WF_mantel)[[i]] <- paste("WF", NAME, sep="")
    }

  WF_mantel_df <- do.call(rbind, WF_mantel)

  print("Smallest Group AVG dist Waterfilter")
  smallest_group_WF <- min(rowSums(otu_table(ps_WF)))
  print(smallest_group_WF)
  avgdist_abund_table_WF  <- avgdist(WF_mantel_df , dmethod = "bray", sample = smallest_group_WF)

  #####################
  #Run the Mantel Test#
  #####################
  dist_tibble_Fish <- avgdist_abund_table_Fish %>%
    as.matrix() %>%
    as_tibble(rownames="sample") %>%
    pivot_longer(-sample)

  dist_tibble_WF <- avgdist_abund_table_WF %>%
    as.matrix() %>%
    as_tibble(rownames="sample") %>%
    pivot_longer(-sample)

  rare_dist_matrix_Fish <- dist_tibble_Fish %>%
    left_join(., metadata, by = c("sample" = "SampleID"))%>%
    left_join(., metadata, by = c("name" = "SampleID")) %>%
    select(sample, name, value) %>%
    pivot_wider(names_from = name, values_from = value) %>%
    select(-sample) %>%
    as.dist()

  rare_dist_matrix_WF <- dist_tibble_WF %>%
    left_join(., metadata, by = c("sample" = "SampleID"))%>%
    left_join(., metadata, by = c("name" = "SampleID")) %>%
    select(sample, name, value) %>%
    pivot_wider(names_from = name, values_from = value) %>%
    select(-sample) %>%
    as.dist()

  abund_temp <-
    mantel(rare_dist_matrix_Fish, rare_dist_matrix_WF, method = "spearman", permutations = 9999, na.rm = TRUE)
  Env.out = data.frame(species=Datasetname, env="WF", mantel=abund_temp$statistic, sig=abund_temp$signif)

    Mantel_list[[Mantel_list_length+1]] <- Env.out
    names(Mantel_list)[[Mantel_list_length+1]] <- paste(Datasetname, "WF", "NoAutumn",sep="_")
}
## [1] "OE"
## [1] "Smalles Group AVG dist"
## [1] 9770
## [1] "Smallest Group AVG dist Waterfilter"
## [1] 3576
## [1] "GC"
## [1] "Smalles Group AVG dist"
## [1] 6441
## [1] "Smallest Group AVG dist Waterfilter"
## [1] 3576
pslist[["Mantel_list"]] <- Mantel_list

# > Mantel_list$`OE_WF_MG-BB`
#   species env    mantel   sig
# 1      OE  WF 0.5307317 1e-04
# > Mantel_list$`GC_WF_MG-BB`
#   species env    mantel    sig
# 1      GC  WF 0.3847558 0.0149

# Mantel_list$`OE_WF_SS-ML`
#   species env    mantel   sig
# 1      OE  WF 0.2997895 1e-04
# Mantel_list$`GC_WF_SS-ML`
#   species env     mantel    sig
# 1      GC  WF 0.07524223 0.1884


#####################
#Plot Mantel Results#
#####################

Mantel_list_clean <- lapply(Mantel_list, function(x) {
  x[complete.cases(x), ]
})
Mantel_out <- as.data.frame(do.call(rbind, Mantel_list_clean))
rownames(Mantel_out) <- NULL
Mantel_out$sig <- as.numeric(Mantel_out$sig)
Mantel_out$mantel <- as.numeric(Mantel_out$mantel)
Mantel_out_sig <- Mantel_out %>% 
  filter(sig < 0.05) 
  
vars  <- c("Salinity", "Temperature", "O2","SPM", "NO3", "NH4", "NO2", "PO4",
"FCF", "GSI", "HSI", "Length", "Age", "WF")   

Envs <- c("Salinity", "Temperature", "O2","SPM", "NO3", "NH4", "NO2", "PO4")
Phy  <- c("Length", "Age", "FCF","HSI", "GSI")

  Mantel_out_sig <- Mantel_out_sig %>% 
    mutate(Species= factor(species, levels = c("OE", "GC", "WF"))) %>%
    mutate(Env=factor(env, levels = vars)) %>%
    mutate(kind = case_when(
      env %in% Envs ~ "Environment",
      env %in% c("WF") ~ "ASV",  
      TRUE ~ "Physiology"))
  Mantel_plot <- Mantel_out_sig %>% 
    ggplot(., aes(x=Env, y=Species, fill=mantel)) +
    geom_tile() +
    scale_fill_gradient2(
      low = "#FFFFFF",
      high = "#006000",
      limit = c(0,1)) +
      theme_minimal() +
      theme(axis.text.x = element_text(angle = 90, size = 15, face ="bold"),
        #axis.title.x.bottom = element_text(size=15,face = "bold"),
        #axis.title.y.left = element_text(size=15,face = "bold"),
        strip.text.x = element_text(size = 15, face = "bold"),
        strip.text.x.bottom = element_text(size = 15,face = "bold"),
        strip.text.y.left = element_text(size = 15,face = "bold"),
        axis.title.x.bottom = element_blank(), 
        axis.title.y.left = element_blank(), 
        axis.text.x.bottom = element_text(face = "bold", angle = 45, hjust = 1),
        axis.text.y.left = element_text(face = "bold", size = 15),
        legend.title = element_text(size = 15, face ="bold"),
        legend.text = element_text(size = 15,face = "bold"),
        plot.title = element_text(size = 15, face = "bold"),
        plot.subtitle = element_text(size = 15),
        plot.caption = element_text(size = 15)) +
    labs(#title = "Correlation relationship between avg.dist Bray-Curtis matrices and varibles", y = "Sample kind", 
      fill="corr") +
      labs(fill = "Mantel's r") +
    #theme(legend.box.background = element_rect(color = "black"))
    theme(legend.position="right", legend.key.width = unit(1, "cm")) 
  Mantel_plot <- Mantel_plot +  geom_text(vjust=0.76,hjust=0.49, aes(label = paste0(c("","*")[(abs(sig) <= .05)+1],
              c("","*")[(abs(sig) <= .001)+1], c("","*")[(abs(sig) <= .001)+1])))
  kindOrder <- c("Environment", "Physiology", "ASV")
  Mantel_plot <- Mantel_plot + facet_grid(. ~ factor(kind, levels = kindOrder), drop = TRUE, scale = "free", space = "free_x")
  
   

  Mantel_plot <- cowplot::plot_grid(Mantel_plot, labels = c(""), ncol = 1)
  Mantel_plot <- plot_grid(Mantel_plot, ncol = 1, rel_heights = c(0.02, 0.05, 0.98))
      plot(Mantel_plot)

ggsave(Mantel_plot, filename = paste(paste("Mantel_test",
          sep="_"),".png", sep=""), path = pathPlots , device='png', dpi=300, width = 7,height = 3)

6.4 Export ASV & Core Dataframe

#############################################
#Create Bacterial Counts and Reseq Dataframe#
#############################################  

SSU_Data_list <- list()
for (Dataset in c(names(pslist)[!grepl("WF", names(pslist))][grepl("ps_GC|ps_OE", names(pslist)[!grepl("WF", names(pslist))])], "ps_WF")) {

  require(plyr)
  require(dplyr)
  require(phyloseq)

  Datasetname <- sub("ps_", "", Dataset)
  ps            <- pslist[[Dataset]]
  SSU_Data_list_length <- length(SSU_Data_list)
  
  
  TAX <- as.data.frame(tax_table(ps))
  OTU <- as.data.frame(t(otu_table(ps)))
  REFSEQ <- refseq(ps)
  REFSEQ <- stack(as.character(REFSEQ, use.names=TRUE))
  colnames(REFSEQ)[colnames(REFSEQ) == "ind"] <- "ASV"

  #################################
  #Creating AverageSums per Season#
  #################################
  SeasonSums <- ps %>%
    transform_sample_counts(function(x) {x/sum(x)*100}) %>%
    phyloseq::otu_table() %>%
    as.data.frame() %>%
    rownames_to_column(var = "SampleID") %>%
    dplyr::left_join(SAMDF16S, by = c("SampleID" = "SampleID")) %>%
    dplyr::group_by(Season) %>%
    dplyr::summarise(dplyr::across(rownames(tax_table(ps)), mean, na.rm = TRUE)) %>% 
    t() %>%
    as.data.frame() %>%
    `colnames<-`(.[1, ]) %>%
    .[-1, ] %>%
    setNames(paste0("avg_", colnames(.))) %>%
    mutate_all(as.numeric) %>%
    rownames_to_column(var="ASV")

    SSU_Data  <- TAX %>%  
     rownames_to_column(var = "ASV") %>% 
       left_join( #Add relative ASVmeans 
       data.frame(t(phyloseq::otu_table(ps %>%
        transform_sample_counts(function(x) {x/sum(x)*100})))) %>%
          dplyr::mutate(ASVmeans = rowMeans(.)) %>%
          dplyr::mutate(ASV = rownames(.)) %>% 
          dplyr::select(ASV, ASVmeans)
       ) %>%
    left_join( #Add Sums by Season
        SeasonSums
       ) %>%
    left_join( #Add relative ASV data per sample
        data.frame(t(phyloseq::otu_table(ps %>%
        transform_sample_counts(function(x) {x/sum(x)*100})))) %>%
        rownames_to_column(var = "ASV")
    ) %>%
     left_join(REFSEQ
    ) %>%
      dplyr::arrange(desc(ASVmeans))
   
   rownames(SSU_Data) <- SSU_Data$ASV

   write.csv2(SSU_Data, paste0(file.path(path_Output_16S, paste(paste(save_name, 
            paste("SSU_ALL", sep=""), Date, Datasetname,  sep="_"),".csv", sep=""))))
   
   ###################  
   #Extract Core-Taxa#
   ###################
    core_taxa_standard <- ps %>%
      transform_sample_counts(function(x) {x/sum(x)*100}) %>%
      microbiome::core_members(., detection = 0.0001, prevalence = 90/100)
    
    #   ps_rel <- microbiome::transform(pslist[[paste("ps_", Species, sep="")]], "compositional")
    # #Initial core computing
    # #Setting a detection threshold of 0.0001% (within samples) and prevelance threshold of 
    #80% (across samples).
    # core_taxa_standard <- microbiome::core_members(ps_rel, detection = 0.0001, prevalence = 90/100)
    
  write.csv2(SSU_Data[SSU_Data$ASV%in% core_taxa_standard,], 
             file.path(path_Output_16S,paste(paste(save_name,"16S_merge_Filter_ASV_besttax_Core", 
                                                   Datasetname, Date, sep="_"), ".csv", sep="")))
  
  SSU_Data_list[[SSU_Data_list_length+1]] <- SSU_Data
  names(SSU_Data_list)[[SSU_Data_list_length+1]] <- Datasetname 
}

#-

Save Data

###########
#Save Data#
###########
saveRDS(pslist, file.path(path_Output_16S, paste(paste(save_name, "ps_16S_merge_Filter_ASV_besttax", Date, sep="_"), ".rds", sep="")))